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These figures show the fate of the product in a transient transport-controlled bimolecular 
reaction under vortex-stirred mixing. The left figure is obtained using a popular numerical 
formulation, which violates the non-negative constraint. The right figure is based on the 
proposed computational framework. These figures clearly illustrate the main contribution of 
this paper: The proposed computational framework produces physically meaningful results for 
advective-diffusive-reactive systems, which is not the case with many popular formulations. 
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species balance for advection-diffusion-reaction equations under the 

finite element method 


M. K. Mudunuru and K. B. Nakshatrala 

Department of Civil and Environmental Engineering, University of Houston. 

Abstract. We present a robust computational framework for advective-diffusive-reactive systems 
that satisfies maximum principles, the non-negative constraint, and element-wise species balance 
property. The proposed methodology is valid on general computational grids, can handle hetero¬ 
geneous anisotropic media, and provides accurate numerical solutions even for very high Peclet 
numbers. The significant contribution of this paper is to incorporate advection (which makes the 
spatial part of the differential operator non-self-adjoint) into the non-negative computational frame¬ 
work, and overcome numerical challenges associated with advection. We employ low-order mixed 
finite element formulations based on least-squares formalism, and enforce explicit constraints on the 
discrete problem to meet the desired properties. The resulting constrained discrete problem belongs 
to convex quadratic programming for which a unique solution exists. Maximum principles and the 
non-negative constraint give rise to bound constraints while element-wise species balance gives rise 
to equality constraints. The resulting convex quadratic programming problems are solved using an 
interior-point algorithm. Several numerical results pertaining to advection-dominated problems are 
presented to illustrate the robustness, convergence, and the overall performance of the proposed 
computational framework. 


1. INTRODUCTION AND MOTIVATION 

Advection-diffusion-reaction (ADR) equations are pervasive in the mathematical modeling of 
several important phenomena in mathematical physics and engineering sciences. Some examples in¬ 
clude degradation/healing of materials under extreme environmental conditions jT], coupled chemo- 
thermo-mechano-diffusion problems in composites [2], contaminant transport [3], turbulent mixing 
in atmospheric sciences j4], diffusion-controlled biochemical reactions [5], tracer modeling in hydro¬ 
geology m, and ionic mobility in biological systems Additionally, ADR equation serves as a 
good mathematical model in numerical analysis, as it offers various unique challenges in obtaining 
stable and accurate numerical solutions [8j. 

The primary variables in these mathematical models are typically the concentration and/or 
the (absolute) temperature. These quantities naturally attain non-negative values. Under some 
popular constitutive models (such as Fourier model and Fickian model, and their generalizations), 
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these physical quantities satisfy diffusion-type equations, which are elliptic/parabolic partial differ¬ 
ential equations (PDEs) and can be non-self-adjoint. These PDEs are known to satisfy important 
mathematical properties like maximum principles and the non-negative constraint (e.g., see [S]). A 
predictive numerical formulation needs to satisfy these mathematical properties and physical laws 
like the (local and global) species balance. It is well-documented in the literature that traditional 
numerical methods perform poorly for advection-dominated ADR equations (e.g., see [SlllOj l. In 
the past few decades, considerable progress has been made to obtain sufficiently accurate numerical 
solutions for ADR equations on coarse computational grids m- It is then natural to ask: “why 
there is a need for yet another numerical formulation for ADR equation?^’. We now outline several 
reasons for such a need. 

(a) Localized phenomena and node-to-node spurious oscillations: For advection-dominated prob¬ 
lems, it is well-known that the standard single-field Galerkin finite element formulation gives 
node-to-node spurious oscillations on coarse computational grids [8]. Moreover, it cannot cap¬ 
ture steep gradients such as interior and boundary layers. Various alternate numerical techniques 
have been proposed to avoid these spurious oscillations |12| . Some methods seem to capture 
steep gradients in interior layers while others capture boundary layers. However, most of them 
do not seem to capture both interior and boundary layers, and at the same time avoid node-to- 
node spurious oscillations m A notable work towards this direction is by Hsieh and Yang m, 
which can capture both interior and boundary layers under adequate mesh refinement. However, 
this formulation has several other deficiencies some of which are discussed below and illustrated 
using numerical examples in subsequent sections. 

(b) Violation of the non-negative constraint and maximum principles: As mentioned earlier, phys¬ 
ical quantities such as concentration and temperature naturally attain non-negative values. It 
is highly desirable for a numerical formulation to respect these physical constraints. This is 
particularly important in a numerical simulation of reactive transport, as a negative value for 
concentration will result in an algorithmic failure. However, it is clearly documented in the 
literature that many existing formulations based on finite element dSHU], finite volume m, 
and finite difference m do not satisfy the non-negative constraint and maximum principles in 
the discrete setting. They also discuss various methodologies to satisfy such properties. But 
most of these methodologies are for pure diffusion equations, which are self-adjoint. For exam¬ 
ple, in Reference |16| . two mixed formulations based on RTO spaces and variational multiscale 
formalism have been modified to meet the non-negative constraint for diffusion equations. This 
approach is not directly applicable to ADR equations for two reasons. First, these formulations 
do not cure the node-to-node spurious oscillations. Second, they do not possess a variational 
structure for ADR equations. Some numerical formulations are constructed to satisfy the non¬ 
negative constraint and maximum principles by taking advection into account (e.g.. 

However, they do not satisfy local and global species balance, and are restricted to isotropic 
diffusion. Conservative post-processing methods exist in the literature to recover certain desired 
properties such as species balance. But such formulations are not variationally consistent |22| . 

(c) Satisfying local and global species balance: In transport problems, the balance of species is an 
important physical law that needs to be met. It is therefore desirable for a numerical formulation 
to satisfy local and global species balance, say, up to machine precision (which is approximately 
10“^® on a 64-bit machine). However, many finite element formulations do not satisfy local 
and global species balance (see [TTlfTill^ f. The main focus of the methods outlined in these 
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papers is to capture the localized phenomena such as boundary and interior layers. Moreover, 
these works did not quantify the errors incurred in satisfying local species balance and global 
species balance. It needs to be emphasized that many finite element methods do exist that 
inherently satisfy local and global species balance, for example, Raviart-Thomas spaces |24) 
and BDM spaces |25| . But these approaches do not fix other issues discussed herein such as 
the node-to-node spurious oscillations or meeting maximum principles for ADR equations. 

(d) Other influential factors: Some other important factors that can affect the performance of 
a numerical formulation are the advection velocity field and its divergence, anisotropy of the 
medium, reaction coefficients, topology of the medium, computational mesh, multiple spatial- 
scales arising due to the heterogeneity of the medium, and multiple time-scales involved in 
various physical processes. Another aspect that brings tremendous numerical challenges is 
chemical reactions involving multiple species. 

It is a herculean task to address all the aforementioned dehciencies, and we strongly believe that 
it may take a series of papers to overcome all the deficiencies. A similar sentiment is shared in the 
review article by Stynes entitled ‘^Numerical methods for convection-diffusion problems or the 30 
years war''' |26| . We therefore take motivation from George Polya’s quote m “ If you can’t solve 
a problem, then there is an easier problem you can solve: find it.’’’ In this paper, we shall pose a 
problem that is simpler than the grand challenge of overcoming all the aforementioned numerical 
deficiencies but still make a significant advancement with respect to the current state-of-the-art. We 
then provide a solution to this simpler problem. To state it more precisely, the main contribution 
of this paper is developing a least-squares-based finite element framework for ADR equations that 
possesses the following properties on general computational grids: 

(PI) No spurious node-to-node oscillations in the entire domain. 

(P2) Captures interior and boundary layers for advection-dominated problems. 

(P3) Satisfies discrete maximum principles and the non-negative constraint. 

(P4) Satisfies local and global species balance. 

(P5) Gives sufficiently accurate solutions even on coarse computational grid^. 

To the best of authors’ knowledge, there exists no finite element methodology for advective-diffusive- 
reactive systems that possesses the desirable properties (P1)-(P5). 

The rest of the paper is organized as follows. Section [2] presents the governing equations for an 
advective-diffusive-reactive system, and discusses the associated mathematical properties. Section 
[3] outlines several plausible approaches, and discusses their drawbacks in meeting the mentioned 
mathematical properties. In Section 01 we propose a constrained optimization-based low-order mixed 
hnite element method to satisfy maximum principles, the non-negative constraint, local species 
balance, and global species balance. In Section [5l we perform a numerical /i-convergence study 
using a benchmark problem. In Section [6l we specialize to transport-limited bimolecular reactions 
to solve problems related to plume development and mixing in isotropic/anisotropic heterogeneous 
media. Finally, conclusions are drawn in Section [71 If one is interested only in the implementation 
of the proposed method, the reader can directly go to Section 01 and appendices A and G. 

We shall denote scalars by lowercase English alphabet or lowercase Greek alphabet (e.g., concen¬ 
tration c and stabilization parameter r). The continuum vectors are denoted by lowercase boldface 
normal letters, and the second-order tensors will be denoted using uppercase boldface normal letters 

^One may expect some subjectivity in calling a mesh to be coarse. But in this paper, we will dehne precisely 
what is meant by a “coarse mesh” for advection-diffusion-reaction equations in terms of M-matrices. 
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(e.g., vector x and second-order tensor D). In the finite element context, the vectors are denoted 
using lowercase boldface italic letters, and the matrices are denoted using uppercase boldface italic 
letters (e.g., vector v and matrix K). We shall use NN to denote non-negative, DMP denotes 
discrete maximum principle, LSB to denote local species balance, and GSB to denote global species 
balance. We shall use XSeed to denote the number of (finite element) nodes in a mesh along x- 
direction. Likewise for YSeed. Other notational conventions adopted in this paper are introduced 
as needed. 


2. GOVERNING EQUATIONS: ADVECTIVE-DIFFUSIVE-REACTIVE SYSTEMS 

Let C be a bounded open domain, where “nd” denotes the number of spatial dimensions. 
The boundary of the domain is denoted by dQ, which is assumed to be piecewise smooth. Mathe¬ 
matically, on := — n, where a superposed bar denotes the set closure. A spatial point is denoted 

by X G n. The gradient and divergence operators with respect to x are, respectively, denoted by 
grad[«] and div[«]. The unit outward normal to the boundary is denoted by n(x). Let c(x) denote 
the concentration field. The boundary is divided into two parts: L'^ and L'^ such that meas(r'^) > 0, 
F U F = dn, and L'^ n L^ = 0. is the part of the boundary on which the concentration is 
prescribed and L'^ is the other part of the boundary on which the diffusive/total flux is prescribed. 

The governing equations for steady-state response of an ADR system take the following form: 


q:(x)c(x) -|- div [c(x)v(x) — D(x)grad[c(x)]] = /(x) 
c(x) = cP(x) 

- Sign[v • n] ^ v(x)c(x) - D(x)grad[c(x)]^ • fi(x) 


in D 
onU 

gP(x) on L'^ 


(2.1a) 

(2.1b) 

(2.1c) 


v(x) is the known advection velocity field, /(x) is the prescribed volumetric source, D(x) is the 
anisotropic diffusivity tensor, q:(x) is the linear reaction coefficient, cP(x) is the prescribed concen¬ 
tration, is the prescribed diffusive/total flux, and the sign function is defined as follows: 


Sign[(^] := 


-1 

if ^9 < 0 


0 

if (/? = 0 

(2.2) 

+1 

if ^9 > 0 



The advection velocity need not be solenoidal in our treatment (i.e., div[v(x)] need not be zero). 
The Neumann boundary condition given in equation (I2.1cl) can be interpreted as follows: 

n(x) • (v(x)c(x) — D (x)grad[c(x)]) = (?*’(x) on L^ (total flux on inflow boundary) (2.3a) 

—n(x) • D(x)grad[c(x)] = (/^(x) on (diffusive flux on outflow boundary) (2.3b) 

where and are, respectively, defined as follows (see Figure [T|): 

r?) := {x G F'^ I v(x) • n(x) < O} (inflow boundary) (2.4a) 

F^ := {x G F'^ I v(x) • n(x) > O} (outflow boundary) (2.4b) 


Remark 2.1. In the literature, more predominantly in the numerical literature, the term ad¬ 
vection is often used synonymously with convection. It should, however, he noted that these two 
terms describe different physical phenomena, and there is a need to clarify the terminology here. An 
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Figure 1. This figure illustrates concentration and flux boundary conditions. Fl corre¬ 
sponds to the inflow boundary while F^ corresponds to the outflow boundary. Total flux is 
prescribed on F^, diffusive flux is prescribed on F^, and concentration is prescribed on F°. 
P = f'^PiF^, Q = F’^PiFI, and R = F^nF^. For well-posedness, we have f'^UF^UFI = dil, 
F'^ n F^ = 0, F^ n F?_ = 0, and F^ n Fl = 0. 


ADR equation arises from the balance of mass of a given species. In ID, an ADR equation takes 
the following form: 


a{x)c{x) -b 


d{vc) 


dx 


A. 

dx 



fix) 


(2.5) 


which is mathematically equivalent to the following equation: 


dv 


dx 


dc 


a(x) + — c{x) + v{x) — 


dx 


jL 

dx 



fix) 


( 2 . 6 ) 


One can obtain a “similar” mathematical equation by linearizing the incompressible Navier-Stokes 
equation, and an appropriate name for this linearized equation is the convection-dissipation-reaction 
(CDR) equation. The CDR equation in ID has the following mathematical form: 


dvQ^ dv 

—v[x) + vq{x) — 
dx dx 


A 

dx 



b{x,po{x)) + 2vo{x) 


dvp 

dx 


(2.7) 


where v{x) is the velocity of the fluid, and po{x) and vo{x) are known pressure and veloeity fields 
about which the Navier-Stokes equation is linearized. From equations and , it is evi¬ 

dent that ID ADR equation and ID CDR equation have similar mathematical forms. However, 
their physical underpinnings are completely different, as the Navier-Stokes equation is obtained by 
substituting a specific constitutive model into the balance of linear momentum. 


2.1. Weak formulations. The following function spaces will be used in the rest of this paper: 

C := {c(x) S H^{D.) I c(x) = c^(x) on P'^} (2.8a) 

W := {tc(x) G H^{n) I t(;(x) = 0 on P'^} (2.8b) 

Q := |q(x) G (L2(r2))"''* | div[q(x)] G T 2 (^^)| (2.8c) 

where q(x) = c(x)v(x) — D(x)grad[c(x)] and H^[Q) is a standard Sobolev space |28| . Given two 
vector fields a(x) and 6(x) on a set K, the standard L 2 inner product over K is denoted as follows: 

{a;b)j^ = j a(x) • 6(x) dif (2-9) 

K 
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The subscript will be dropped '\i K = Q,. The most popular way to construct a weak formulation is to 
employ the Galerkin formalism. Based on the manner in which one applies the divergence theorem, 
the single-field Galerkin formulation for equations (|2.1al) - (l2.1cp can be posed in two different ways. 
Single-field Galerkin formulation 7^1 (SGi): Find c(x) G C such that we have 


{w; ac) — (grad[rc] • v; c) -|- (grad[rt;]; D(x)grad[c]) -|- ( w, 


1 + Sign[v • n] 


n c 


r9 


= {w;f) - V'w;(x) G W (2.10) 

Single-field Galerkin formulation ^2 (SG 2 ): Find c(x) G C such that we have 

1 — Sign[v • n] 


(re; (a-|-div[v]) c) -1- (tc;grad[c] • v) -|- (grad[u;]; D(x)grad[c]) — I w; 


(v • n) c 


r9 


= (w; f) - (w; qP)r, Vw(x) G W 


( 2 . 11 ) 


Note that the solution obtained will be the same regardless whether we use either SGi or SG 2 . 
However, this is not true if one uses total/diffusive flux on Neumann boundary without giving 
due consideration to inflow and/or outflow Neumann boundary conditions. For more details, see 
subsection 12.31 


2.2. Maximum principles and the non-negative constraint. A basic qualitative property 
of elliptic boundary value problems is the maximum principle. This property gives a priori estimate 
for c(x) in fi through its values on F'^. The following assumptions will be made to present a 
continuous maximum principle for ADR equations with both Dirichlet and Neumann boundary 
conditions: 

(Al) D is piecewise smooth domain with Lipschitz continuous boundary dQ. 

(A2) The scalar functions a : D —)• R, (v)j : D —^ R, and (D)jj : D —)• R are continuously 
differentiable in their respective domains for i = 1, • • • ,nd. Furthermore, / G L 2 (D), q'^ G 
L 2 (r'?), and = g*\T'^ with g* G Ff^(D). 

(A3) The diffusivity tensor is assumed to be symmetric, uniformly elliptic, and bounded above. 
That is, there exists two constants (i.e., independent of x), 0 < 71 b < 7 ub < + 00 , such that 
we have 

0 < 7 iby y < y D(x)y < 7 ,by y Vy g R’"^{o} ( 2 . 12 ) 

(A4) The advection velocity field v(x) and the reaction coefficient a(x) are restricted as follows: 

0 < a(x) -|- div [v(x)] < /3i(x) Vx G D (2.13a) 

0 < a{x) -|- 2 *^^^ [''^(^)] ^ / 52 (x) V X G D (2.13b) 

0 < |v(x) • S(x)| < Psix) Vx G n (2.13c) 

where /3i(x) G Lnd/ 2 i^), ; 02 (x) G Lnd/ 2 {^), and Ps{x) G L„d-i(r‘'). It is assumed that 
functions /3i(x), / 32 (x), and Psix) are bounded above for a unique weak solution to exist 
based on the Lax-Milgram theorem. 

(A5) The part of the boundary on which Dirichlet boundary condition is enforced is not empty 
(i.e., T- / 0 ). 

We shall use the standard abbreviation of a.e. for almost everywhere |28j . 
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Theorem 2.2 (A continuous maximum principle). Let assumptions (A1)-(A5) hold and let 
the unique weak solution c(x) of equations (|2.1ap - (I2.1cp belong to C'^(ri) H (7*^(0). If /(x) £ L2{0,) 
and gP(x) G L 2 (r'^) satisfy: 

/(x) < 0 a.e. in Ll (2.14a) 

f 7 ’’(x) > 0 a.e. on r!() U (2.14b) 

then c(x) satisfies a continuous maximum principle of the following form: 


[c(x)] < max 


0, max ic'’(x)l 
xGr= ^ ■' 


In particular, if c^{x.) > 0 then 


m^ [c(x)] = max [cP(x)] 
xgh xGr= 


//cP(x) < 0 then we have the following non-positive property: 

m^ [c(x)] < 0 
xGfl 

Proof. Let and m(x) are, respectively, defined as follows: 

‘f’max := max 0, max [cP(x)] 

xGF'^ 

m(x) := max [0, c(x) - 4>max] 


(2.15) 

(2.16) 

(2.17) 

(2.18) 
(2.19) 


It is easy to check that m(x) is a non-negative, continuous, and piecewise (7^(11) function. From 
equation (|2.19p . it is evident that m(x)|rc = 0 and c(x) = m(x) -|- ^max for any x G 14 unless 
m(x) = 0. Moreover, m(x) G W. By choosing w(x) = m(x), the weak formulation given by 
equation (|2.10p becomes: 


(m; a(m -h 4'max)) - (grad[m] • v; (m 4>max)) + (grad[m]; D grad[m]) 

' 1 -|- Sign[v • fi] 


+ m: 


V • n (m -h 


F'j 


= (m;/) - (m;gP)r9 (2.20) 

It is easy to establish the following identities: 

(m;vn(m-l- 4>max))r-j = (grad[m] • v; (m -P 4>max)) + (m; div[v] (m -h 4>max)) 

-|- (m; grad[m] • v) (2.21a) 

2(grad[m] • v; (m-h 4>max)) = (m;vn(m-h 4>max))r9 - (m; div[v] (m -j- 4>max)) 

- (^max; grad[m] • v) (2.21b) 

(^maxi grad[m] • v) = (^maxi V • n m)r-j - ($max; div[v] m) (2.21c) 


(grad[m] • v; (m -P 4>max)) = ( m; v • n ( 4>max + 


r? 


- { m; div[v] { 4>max + -jm 


(2.21d) 
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Using the above identities, equation (12.201) can be written as follows: 


m; [ a + -div[v] ) m ) + (m; (a + div[v]) $max) + (grad[m]; D grad[m]) 


V • n 

+ ( m;-ml — | m; 

ri 


1 — Sign[v • n] 


^ (v • n) d>max^ = (m-, f) - (m; gP)p, (2.22) 


From equations (12.121) and (I2.13ah - (l2.13cj) . it is evident that 


m; ( a + -div[v] ) m ) + (m; {a + div[v]) <hmax) + (grad[m]; D grad[m]) 


V • n 

+ ( "i; — 2 — j - ( m; 


1 — Sign[v • n] 


V • n 


> 0 


ri 


From equation (I2.14D we have: 


Therefore, one can conclude that 


1 


(m;/) - (m;gP)r, < 0 


(2.23) 

(2.24) 


m; ( a + -div[v] ) m ) + (m; (a + div[v]) d'max) + (grad[m]; D grad[m]) 


V • n 

+ I m; —-— m 




— m 


1 - Sign[v • n] \ _ \ 

-i (v • n) $max 


= 0 


/ T'J 


(2.25) 


In the light of assumption (A3) and equation (12.251) . we need to have grad[m] = 0 (as D(x) is 
bounded below by a constant 7ib)- This further implies the following: 


m(x) = (/>o > 0 Vx G n 


(2.26) 


where <j)Q is a non-negative constant. Since m(x)|rc = 0 and meas(r^) > 0, we have (j)Q = 0. This 
implies that c(x) < 4>inax) which further implies the validity of the inequality given by equation 
(j2.15l) . Finally, equations (I2.16P and (|2.17p are trivial consequences of equation (12.151) . □ 

We have employed the SGi formulation in the proof of Theorem 12.21 One will come to the same 
conclusion even under the SG 2 formulation. By reversing the signs in equation (12.141) . one can easily 
obtain the following continuous minimum principle. 

Corollary 2.3 (A continuous minimum principle). Let assumptions (A1)-(A5) hold and let 
the unique weak solution c(x) of equations (12.lap - (I2.1cp belong to C'^(n) H C®(U). If /(x) G L 2 {^) 
and (?P(x) G L 2 (F'^) satisfy 


/(x) > 0 a.e. in U 
(/'’(x) < 0 a.e. on F^ U F^ 

then c(x) satisfies a continuous minimum principle of the following form: 


min [c(x)] > min 
xeo 


0, min [cP(x)l 
xer^: ^ 


(2.27a) 

(2.27b) 

(2.28) 


In particular, if c^{x.) < 0 then 


min[c(x)] = min [cP(x)] 


(2.29) 




































//cP(x) > 0 then we have the following non-negative property: 

min [c(x)] > 0 (2.30) 

This paper also deals with transient analysis, and the details are provided in Sections |4] and [6l 


2.3. On appropriate Neumann BCs. Many existing numerical formulations m and pack¬ 
ages such as ABAQUS |30| . ANSYS |31| . COMSOL |32j . and MATLAB's PDE Toolbox |33| do not 
pose the Neumann BCs in correct form for advection-diffusion equations. These formulations and 
packages either use the diffusive flux or the total flux on the entire Neumann boundary without 
discerning the following situations: 

• Do we have inflow (i.e., v • n < 0) on the entire Neumann boundary? 

• Do we have outflow (i.e., v • n > 0) on the entire Neumann boundary? 

• Or do we have both inflow and outflow on the Neumann boundary? 

These conditions will dictate whether the resulting boundary value problem is well-posed or not. 
If a numerical formulation does not take into account these conditions, the numerical solution can 
exhibit instabilities, which will have dire consequences in mixing problems. To illustrate, consider 
the following ID boundary value problem: 

-^(vc — D^]=0 'ix£(0,L) (2.31a) 

ax \ ax J 

c{x = 0) = Co (2.31b) 


where v, D and cq are constants, and L is the length of the domain. We now consider two different 
cases for the Neumann BC: 

^uc — = Qo Sit x = L (2.32a) 

dc 

—D— = qo Sit x = L (2.32b) 

ax 

where qq is a constant. Equation (j2.32aj) corresponds to the total flux BC while equation ()2.32b(l is 
the diffusive flux BC. The analytical solutions for these two different Neumann BCs are, respectively, 
given as follows: 


Ci{x) = ^ (^go + {vco - go) e o ^ 

1 / -vL v{x-L) \ 

C 2 {x) = - (^uco + qoe D - q^e d j 


(2.33a) 

(2.33b) 


The solution ci(x) blows if u > 0, and C 2 (x) blows if u < 0. On the other hand, the exact solution 
based on the Neumann BC given in equation (I2.1cp is well-posed for both inflow and outflow cases. 

To summarize, the boundary value problem is well-posed under the prescribed diffusive flux on 
the entire Neumann boundary if the flow is outflow on the entire r*?. The boundary value problem 
is well-posed under the prescribed total flux on the entire Neumann boundary if the flow is inflow on 
the entire T*?. The Neumann BC given by equation (|2.1cl) is more general, and the boundary value 
problem under this BC is well-posed even if the Neumann boundary is composed of both inflow and 
outflow. 
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3. PLAUSIBLE APPROACHES AND THEIR SHORTCOMINGS 


There are numerous numerical formulations available in the literature for advective-diffusive- 
reactive systems. A cavalier look at these formulations can be deceptive, as one may expect more 
than what these formulations can actually provide. We now discuss some approaches that seem plau¬ 
sible to satisfy the maximum principle and the non-negative constraint for an advective-diffusive- 
reactive system, and illustrate their shortcomings. This discussion is helpful in two ways. First, 
it sheds light on the complexity of the problem, and can bring out the main contributions made 
in this paper. Second, the discussion can provide a rationale behind the approach taken in this 
paper in order to develop the proposed computational framework. To start with, it is well-known 
that the single-field Galerkin formulation does not perform well, as it produces spurious node-to- 
node oscillations on coarse grids m The formulation also violates the non-negative constraint 
and maximum principles for anisotropic medium, and does not possess element-wise species balance 
property [MITT]. 

3.1. Approach Clipping/cut-off methods. There are various post-processing proce¬ 
dures such as clipping/cut-off methods [2^1^ to ensure that a certain numerical formulation 
satisfies the non-negative constraint. The key idea of these methods is to simply chop-off the nega¬ 
tive values in a numerical solution. Although a clipping method is a variational crime, this approach 
appeals the practitioners because of its simplicity. However, there are many reasons, which are often 
overlooked by the practitioners, why a clipping method is not appropriate for ADR equations with 
anisotropic diffusivity. The reasons, which are documented below, go well beyond the philosophical 
issue of “variational crime.” The reasons should sufficiently justify and motivate to employ a rather 
sophisticated computational framework just like the one proposed in this paper. 

(i) The violation of the non-negative constraint is small only for pure diffusion equations with 
isotropic diffusivity. The violations can be large in the case of anisotropic diffusion. If the 
maximum eigenvalue is not much smaller than unity, then naive /i/p-refinement will not always 
reduce the negative values and clipping procedure can give erroneous results. Figure 19 and 
problem 6.2 in the paper illustrate this point. This has been illustrated even for diffusion 
equations in Reference |35j . 

(ii) Although tensorial dispersion frequently arises in the modeling of subsurface systems, many 
practitioners employ isotropic diffusion in their numerical simulations just to avoid large non¬ 
negative violations in their reactive-transport modeling. As mentioned earlier, in the case 
of isotropic diffusion, one can go away with the clipping procedure. But there is a need 
for predictive simulations for realistic scenarios (e.g., anisotropic diffusivity), and one needs 
carefully designed computational frameworks. Simple approaches like the clipping procedure 
will not suffice. 

(iii) A clipping procedure, by itself, does not ensure local species balance. 

(iv) The clipping procedure cannot eliminate the spurious node-to-node oscillations. 

(v) The ramifications of clipping the negative values on the species balance and on the overall 
accuracy of solutions have not been carefully studies or documented. 

(vi) Finally, both h- and p-refinements may decrease the negative values and reduce spurious node- 
to-node oscillations for advection-dominated and reaction-dominated ADR problems. How¬ 
ever, our objective is to satisfy maximum principles, non-negative constraint, species balance, 
reduce spurious node-to-node oscillations, and obtain sufficiently accurate numerical solutions 
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on coarse computational grids. Extensive mesh and polynomial refinements defeats the main 
purpose, as these approaches will incur excessive computational cost. 


3.2. Approach ^2: Mesh restrictions. Recently, there has been a surge on the study of 
constructing meshes to satisfy various discrete maximum principles both within the context of 
single-field and mixed finite element formulations |36H38j . The primary objective of these methods 
is to develop restrictions on the computational meshes to meet the underlying principles. However, 
it should be noted that there are various drawbacks for these methods. The important ones are 
described as follows: 

(i) Most of these mesh restriction methods are for simplicial meshes (such as three-node triangular 
element and four-node tetrahedral element). Extending these results to non-simplicial elements 
is not trivial or may not be possible. 

(ii) The boundary conditions are restricted to only Dirichlet on the entire boundary of the domain. 
Incorporating mixed boundary conditions or a general Neumann BC given by equation (I2.1cp 
has not been addressed. 

(iii) Generating a DMP-based mesh for complex domains is extremely difficult and sometimes 
impossible. 

(iv) For highly advection-dominated and reaction-dominated problems, we need a highly refined 
DMP-based meshes. Constructing such meshes is computationally intensive. 

(v) Even though the mesh restriction conditions put forth for the weak Galerkin method by Huang 
and Wang m is locally conservative, it is restricted to pure anisotropic diffusion equations. 
Generalizing it to obtain locally conservative DMP-based meshes for anisotropic ADR equa¬ 
tions is not apparent. Moreover, it still suffers from the above set of drawbacks. 


3.3. Approach #3: Using non-negative methodologies for diffusion equations. Re¬ 
cently, optimization-based finite element methods are proposed to satisfy the non¬ 

negative constraint and maximum principles for diffusion-type equations. These non-negative 
methodologies are for self-adjoint operators and are constructed by invoking Vainberg’s theorem |39| . 
That is, they utilize the fact that there exists a scalar functional such that the Gateaux variation 
of this functional provides the weak formulation and the Euler-Lagrange equations provide the 
corresponding governing equations for the diffusion problem. Corresponding to this continuous 
variational/minimization functional, a discrete non-negative constrained optimization-based finite 
element method is developed. Unfortunately, such a variational principle based on Vainberg’s theo¬ 
rem does not exist for the Galerkin weak formulation for an ADR equation, as the spatial operator 
is non-self-adjoint |40| . 


3.4. Approach Posing the discrete equations as a P-LCP. Let h be the maximum 
element size, ||v||cxD,r 2 be the maximum value for advection velocity field, Ooo.o be the maximum 
value for linear reaction coefficient, and Amin be the minimum eigenvalue of D(x) in the entire 
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domain. Mathematically, these quantities are defined as follows: 


h := 


max 


max [|(v(x))i 

l<i<nd 


I V||oo,f2 

Q:cx),r2 

^min ■= [^min,D(x)] 

A max 


Vx G n 


max [q:(x)] 

xGH 


[^max,D{x)] 


(3.1a) 

(3.1b) 

(3.1c) 

(3.1d) 

(3.1e) 


where Qh is a regular linear finite element partition of the domain Q such that Qh = {}e=i^e- “A^e/e” 
is the total number of discrete non-overlapping open sub-domains. The boundary of fie is denoted 
as clfle := fie ~ fie- is the diameter of element fig. Amin.D(x) Ainax,D(x) respectively, the 
minimum and maximum eigenvalue of D(x) at a given point x G fl. Correspondingly, the element 
Peclet number Pe^ and the element Damkohler number Da/i are defined as follows: 

Pe/, := (3.2a) 

Da,, := (3.2b) 

^min 

Herein, Da,, is defined based on linear reaction coefficient and diffusivity. However, it should be 
noted that there are various ways to construct different types of element Damkohler numbers (for 
instance, see Reference |41) for isotropic diffusivity). 

After low-order finite element discretization of either SGi or SG 2 , the discrete equations for the 
ADR boundary value problem take the following form: 


Kc = f (3.3) 

where K is the stiffness matrix (which is neither symmetric nor positive definite), c is the vector 
containing nodal concentrations, and f is the volumetric source vector. The matrix K is of size 
ncdofs X ncdofs, where “ncdofs” denotes the number of free degrees-of-freedom for the concentra¬ 
tion. The vectors c and f are of size ncdofs x 1. 

In the rest of this paper, the symbols ^ and ^ will be used to denote the component-wise 
comparison of vectors and matrices. That is, given two vectors a and b, a ^ b means that (a )i < {b)i 
for all i. Likewise, given two matrices A and B, A < B means that {A)ij < {B)ij for all i and j. 
The mathematical means of the symbols -< and should now be obvious. We shall use 0 and 
O to denote zero vector and zero matrix, respectively. 


Definition 3.1 (P-matrix, Z-matrix, and M-matrix). A matrix A G is a P- 

matrix if ^ (A-|-A"*") is positive-definite. The matrix is a Z-matrix if (A)ij < 0, where i j 
and i, j = 1, ■ ■ ■ ,nd. The matrix is an M-matrix if A is a P-matrix and a Z-matrix. 


Definition 3.2 (Coarse mesh demarcation for anisotropic ADR equations). A regular low- 
order finite element computational mesh is said to he eoarse with respect to 
(a) spurious oscillations f/Pe^ > 1 

(h) spurious oscillations and large linear reaction coefficient i/Pe,, > 1 and Da,, > 1 
(c) spurious oscillations, large linear reaction coefficient, and a discrete maximum principle if the 
stiffness matrix K associated with either SGi or SG 2 is not an M-matrix 
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It can be easily shown through counterexamples that the stiffness matrix K for ADR equation 
will not always be a Z-matrix. We shall now provide two such counterexamples. The first coun¬ 
terexample is the low-order finite element discretization based on two-node linear element for the 
following ID ADR equation (with constant velocity, diffusivity, and linear reaction coefficients): 

dc cl‘^ c 

ac + V- - D——^ = fix) Vx € D := (0,1) (3.4a) 

ax dx^ 

c{x) = c^{x) Vx e (9D := {0,1} (3.4b) 


with a > 0, D > 0, and u € M. The entries of stiffness matrix K for an intermediate node (using 
equal-sized two-node linear finite element) is given as follows: 


ah 

~6 


[1 4 1] 


11 

[ -r 1 J 

f C-l ] 

I P r 1 J 

f 1 

Ci 


1 1 



1 Ci+1 J 

1 1 

1 Ci +1 ) 

1 Ci +1 ) 


(3.5) 


On trivial manipulations on equation (13.5p . it is evident that the stiffness matrix is a Z-matrix if 
and only if the following condition is satished: 


h <h 


max 


12D 

3|u| -|- \/9u^ + 24aD 


(3.6) 


which is not always the case. The second counterexample is based on a simplicial finite element 
discretization (e.g., three-node triangular/four-node tetrahedron element) of ADR equation with 
Dirichlet BCs on the entire boundary. If any nd-simplicial mesh does not satisfy the following 
condition then K is not a Z-matrix |38l Theorem 4.3]: 


0 < 


hn 


+ 




< cos(/3 ~-i) 


{nd + 1) {nd + 1) {nd + 2) 

Vp, g = 1, 2, • • • ,nd+l, pj^q, ^ 


(3.7) 


where p and q are the any two given arbitrary vertices of Dg. Due is the integral element average 
anisotropic diffusivity. A^.^ denotes the minimum eigenvalue of Dn^ ■ hp and hg are the 
perpendicular distance (or the height) from vertices p and q to their respective opposite faces. 
/3 ^-1 is the dihedral angle measured in metric between two faces opposite to vertices p and 
q of an element Dg. 


Proposition 3.3 (P-matrix linear complementarity problem for ADR equations). Given 
that assumptions (Al)-(A5) hold, then the stiffness matrix K associated with low-order finite 
element discretization of either SGi or SG 2 is a P-matrix. Furthermore, if c has to be DMP- 
preserving on any arbitrary coarse mesh, then the resulting constrained discrete equations of single¬ 
field Galerkin formulation can be posed as a P-LCP. 


Proof. We prove only for SGi formulation and extending it to SG 2 is a trivial manipulation. 
The symmetric part of the element stiffness matrix is given as follows: 


1 

2 


{Ke + Kj) 


= j (a(x) + ^div[v(x)]^ 

Oe 


N^N dDe + / PD(x)P^ dDe 

He 

-h j dP^ 

rl 


(3.8) 
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where N is row vector containing shape functions and B = {DN)J ^ (see Appendix A for details 
on DN and J). From equation (|3.8p and assumptions (A1)-(A5), it is evident ^ [Kf, + KJ) is 
positive semi-definite. Furthermore, the assembly procedure ensures that the global stiffness matrix 
K is positive definite |42l Section 2 and Section 3]. As the mesh is coarse, K is not an M-matrix. 
But we want c to satisfy the DMP constraints. Hence, this results in the following constrained 
discrete system of equations: 

Kc = f + X (3.9a) 

A ^ 0 (3.9b) 

c ^ 0 (3.9c) 

A*c = 0 (3.9d) 

As FC is a P-matrix, the above system is a P-matrix linear complementarity problem. This completes 
the proof. □ 

It needs to be emphasized that solving a LCP problem with P-matrix is, in general, NP-hard |43| . 
Therefore, posing the discrete equations as a LCP problem and numerically solving it is not a 
viable approach, especially for large-scale ADR problems with high Pe/j. Moreover, it is not always 
feasible to find a DMP-based /i-refined mesh that will produce accurate results for ADR equation 
for sufficiently high element Peclet number and element Damkbhler number. 


3.5. Approach ^5: Posing the discrete problem as constrained normal equations. 

One way of constructing an optimization-based approach to meet the non-negative constraint is to 
rewrite the discrete problem as the following constrained normal equations: 


minimize - {c\ Kc) — {c] f) 

C e R-r^t^d-ofs 2 

(3.10a) 

subject to c P 0 

(3.10b) 

where (•;•) denotes the standard inner-product in Euclidean spaces, 
optimality conditions can be written as: 

The corresponding first-order 

K^Kc = K^f + A 

(3.11a) 

c P 0 

(3.11b) 

A P 0 

(3.11c) 

XiCi = 0 Vi = 1, • • • , ncdofs 

(3.11d) 

If there are no constraints, the optimization problem becomes: 


minimize - {c, Kc) — {c; f) 

c e 2 

(3.12) 

The first-order optimality condition for the unconstrained discrete optimization problem is: 

K^Kc = K^f 

(3.13) 

In the numerical mathematics literature fe.g.. see |44jL the above svstem of equations (]3.13p is 


referred to as normal equations. The three main deficiencies of this approach are: 

(i) The constrained optimization-based normal equations method does not avoid node-to-node 
spurious oscillations. In addition, there is no obvious way of fixing the method to avoid this 
type of unphysical solutions. 
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(ii) It is well-known that the condition number of K^K will be much worse than K. So the 
numerical solution will be less reliable, less accurate, and numerically not stable |44| . 

(hi) The discrete optimization problem given by equation (13.121) on which non-negative constraints 
are enforced does not have a corresponding continuous variational/minimization problem. 

We shall use the following academic problem to illustrate the aforementioned deficiencies: 

'/x£(0,l) (3.14a) 

c{x = 0) = c{x = 1) = 0 (3.14b) 

where v, D, and / are assumed to be constants. In our numerical experiment, we have taken the 
number of mesh elements to be 11, v/D = 150, and / = 1. The element Peclet number (Pe^ = ^) 
is approximately 6.82. Since it is greater than unity, there be will spurious node-to-node oscillations 
under the standard single-field Galerkin formulation. From Figure [2l it is evident that the normal 
equations approach does not eliminate the spurious node-to-node oscillations. The condition number 
of the stiffness matrix under the standard single-field Galerkin formulation is 8.41, whereas the 
condition number of the stiffness matrix under the normal equations approach is 70.69. For small 
element Peclet numbers, deficiency (i) can be avoided. But dehciencies (ii) and (iii) will still be 
present and cannot be circumvented. Hence, posing the discrete problem as constrained normal 
equations is not a viable approach to meet maximum principles and the non-negative constraint. 

4. PROPOSED COMPUTATIONAL FRAMEWORK 

We employ least-squares formalism to develop a class of structure-preserving numerical formula¬ 
tions whose solutions satisfy DMP, LSB, and GSB. The formulations are built based on minimization 
of unconstrained/constrained quadratic least-squares functionals. In a least-squares-based finite el¬ 
ement formulation, a non-physical least-squares functional is constructed in terms of the sum of 
the squares of the residuals in an appropriate norm. These residuals are based on the underlying 
governing equations. However, it should be noted that LSFEMs are different from the Galerkin 
least-squares or stabilized mixed methods, where least-squares terms are added locally or globally 
to variational problems. 

The success of LSFEM is due to the rich mathematical foundations that influence both the anal¬ 
ysis and the algorithmic development. LSFEM offers several attractive features. The resulting weak 
formulations are coercive. Hence, a unique global minimizer exists for the least-squares functional 
and this minimizer coincides with the exact solution. Gonforming finite element discretizations of 
least-squares functionals leads to stable and (eventually) optimally accurate numerical solutions. 
For mixed LSFEM-based formulations, equal order interpolation can be used for all the unknowns, 
which is computationally the most convenient. The resulting algebraic system is symmetric and pos¬ 
itive definite. Thus, the discrete system can be solved using standard and robust iterative numerical 
methods. For more details on LSFEM for various applications, see Bochev and Gunzberger |45| 
and Jiang |46| . 

4.1. Design synopsis of the proposed numerical methodology. The central idea of the 
proposed computational framework is to constrain a least-squares functional with LSB and non¬ 
negative constraints. The main steps involved in the design of the proposed computational frame¬ 
work are: 

(i) The governing equations of the ADR problem are written in first-order mixed form. 
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(ii) We construct a stabilized least-squares functional for these first-order governing equations. 

(iii) We construct algebraic equality constraints to enforce element-wise/local species balance (LSB). 

(iv) We enforce bound constraints to the constructed LSFEM to meet maximum principles and the 
non-negative constraint in the discrete setting. In order to achieve this, we shall use low-order 


finite element interpolation for c(x). 

The first-order mixed form of the governing equations can be written as: 

q(x) — v(x)c(x) -|- D(x)grad[c(x)] =0 in (4.1a) 

div[q(x)] = /(x) — q;(x)c(x) in (4.1b) 

c(x) = cP(x) on (4.1c) 

(^q(x) - + v(x)c(x)^ • a(x) = gP(x) on n (4.1d) 

The bound constraints to meet discrete maximum principles take the following form: 

Cminl C ^ Cmaxl in (4'2) 


where Cmin and Cmax are the minimum and maximum concentration values possible in 11. The LSB 
equality constraints can be constructed in two different ways. The first approach is based on the 
integral statement of the balance of species on an element, and takes the following mathematical 
form: 

J a(x)c(x) dlle + J q(x) • n(x) dTe = J /(x) dflg (4.3) 

The second approach is to enforce equation (I4.1bp in each mesh element fie in an integral form: 

J a(x)c(x) dOe -h J div[q(x)] dfle = J /(x) dflg (4.4) 

One can obtain equation (|4.3I) by applying the divergence theorem to equation (|4.4p . which means 
that these two approaches are equivalent in the continuous setting. This will not always be the 
case in the discrete setting. In the case of simplicial and non-simplicial low-order finite elements, 
these approaches are equivalent. However, these two approaches can be different in the case of 
higher-order finite elements. This is because in certain higher-order finite elements (e.g., nine-node 
quadrilateral element), not all the nodes are on the boundary of the element. The flux at an interior 
node contributes to the second integral in equation (j4.4p but not to the corresponding term in 
equation (14.3p . These issues are beyond the scope of this paper. Herein, we take the first approach 
given by equation (14.3p . 

We next construct two different least-squares functionals and analyze the influence of various 
constraints on the performance of these LSFEMs. It should be noted that Hsieh and Yang m have 
proposed similar least-squares functionals, but they considered homogeneous isotropic steady-state 
advection-diffusion equations. Moreover, even in the simple setting of isotropic diffusivity, they did 
not consider general Neumann BCs, spatially varying velocity fields, simplicial vs. non-simplicial 
elements, or the effects of DMPs and LSB on the performance of the least-squares functionals. This 
paper investigates all the mentioned aspects: we incorporate anisotropy, heterogeneity, transient 
effects, linear reaction terms, non-solenoidal spatially varying velocity fields, and DMP and LSB 
constraints. 
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4.2. Weighted primitive LSFEM. The weighted primitive LSFEM is the standard way of 
constructing a LSFEM-based formulation. It does not contain any additional stabilization terms. 
The weighted primitive least-squares functional 5prim(c, q) ; C x Q —)• M in L 2 -norm can be written 
as: 


1 2 
S'Prim (c, q) := - A(x) (q - cv -h Dgrad[c]) 




1 


/3(x)(ac + div[q] - /) 




/1 + Sign[v • n] \ 

V 2 J 



• h — qP 


2 

r-j 


(4.5) 


where the second-order tensor A(x) and the scalar function /3(x) are the weights, which are defined 
as follows: 


A(x) = j 

[ D-V2(x) 

LS Type-1 

LS Type-2 


(4.6a) 

|5(x) = 1 


1 

if a(x) = 01 
if a(x) / OJ 

LS Type-1 

LS Type-2 

(4.6b) 


A corresponding weak form can be obtained by setting the Gateaux variation of the functional (14.5p 
to zero. We shall show in Sections[5]and|6]that a naive way of constructing LSFEM formulation, just 
like the weighted primitive LSFEM, does not perform well for advection-dominated ADR problems. 
Moreover, enforcing LSB and DMP constraints do not seem to have a profound effect. In order 
to adequately capture steep boundary and interior layers, we introduce an alternate stabilized 
LSFEM formulation, which will be referred to as the weighted negatively stabilized streamline 
diffusion LSFEM. This stabilized LSFEM formulation will be able to handle a wide spectrum of 
ADR problems ranging from advection-dominated to reaction-dominated problems. 


4.3. Weighted negatively stabilized streamline diffusion LSFEM. The underlying idea 
of the proposed stabilized LSFEM formulation is to combine the streamline diffusion and stabilized 
Galerkin formulations. This is motivated by the prior studies that combining these two formulations 
exhibit enhanced stability (for example, see |14(I47| ). In this formulation, we introduce a small 
element-wise stabilization parameter to correct q(x) in the streamline direction by adding 
second-order derivatives of c(x). The modihed flux along the streamline direction takes the following 
form: 

q = cv — Dgrad[c] -|- (div[cv — Dgrad[c]]) (4-7) 

The basic philosophy of the correction to the flux given by equation (14.7p is in the spirit of stabilized 
finite element formulations such as SUPG |47[I48| . This flux correction is different from that of 
the Flux-Corrected Transport (FCT) methods |49| . Correspondingly, the species balance equation 
accounting for these corrections will be: 

ac + divlq] = f + (4.8) 

where 

/fee •= (grad[/ - ac] • V -b div[v] (/ - ac )) 
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(4.9) 















The modification to the flux (given by equations (I4.7p - (l4.9l) i will present two different ways of 
constructing Neumann BCs. 

The first way utilizes the quantities q(x), c(x), a(x), and /(x), and takes the following form: 

1 + Sign[v —^ (' 4 gQ) 



The second way utilizes q(x), c(x), and the first and second derivatives of c(x). The corresponding 
expression for Neumann BCs takes the following form: 

^ cv — (div[cv — Dgrad[c]]) • n(x) = q^{x.) on T'^ (4.11) 



In the continuous setting, equations (|4.10p and (14.111) are equivalent. However, in the discrete set¬ 
ting, the performance of these equations can be different based on the kind of (finite) element being 
employed. For example, for simplicial elements (such as three-node triangular (T3) element and 
four-node tetrahedral (T4) element) and four-node quadrilateral (Q4) element, both div[grad[c(x)]] 
and grad[grad[c(x)]] are zero for He G T'^. This is because the Hessian of N, which is DDN, is 
a zero matrix for both two-node linear (L2) and three-node triangular elements. For more details, 
see Appendix A. But, this is not the case with non-simplicial linear finite elements for nd = 3 and 
higher-order finite elements (in any dimension). Hence, the Neumann BCs based on equation (14.lip 
are not always physically consistent. However, Neumann BCs based on equation (|4.10p are always 
consistent irrespective of the finite element used. Herein, we have chosen Neumann BCs given by 
equation (I4.10p . 

Based on the above set of equations (I4.7p - (|4.10p . we construct a L 2 -norm based least-squares 
functional. Additionally, as in the Galerkin least-squares method, we add a stabilization term to 
this functional. This stabilization term is as follows: 


1 

2 


E 




div [cv — 


Dgrad[c]] -|- ac — / 


2 


fie 


(4.12) 


The least-squares functional for the weighted negatively stabilized streamline diffusion formulation 
i?NgStb(C) q) : C X Q — )• M in L 2 -norm takes the following form: 


1 


5'NgStb (c, q) := g ^ A(x)(q - cv + Dgrad[c] - 6q^v (div[cv - Dgrad[c]])) 










fiesr^ 


q 


1 -|- Sign[v • n] 


cv -6n^{f -ac)v] • n - q 


+ X ^ div[cv - Dgrad[c]] -k ac - / 




(4.13) 
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The element dependent parameters < 0 and < 0 are given as: 

= -- 




^Lx + <^1 max 

xgfl 


(a + div[v])^ /i^ + S 2 m^ [||div[D] |p] 
-I xGH 


= - 


T \2 u2 


(a + div[v])^ /i^ + r 2 m^ [||div[D]||^l h" 

-I xefi 


^max + n max 


(4.14a) 


(4.14b) 


where 60 , 5i, S 2 , To, Ti, and T 2 are non-negative constants. Appendix B provides a thorough 
mathematical justification behind the above stabilization parameters. 

For unconstrained LSFEMs, the errors incurred in satisfying LSB and GSB can be calculated 
as: 


Je) _ 

^LSB — 


J q;(x)c(x) dn -|- J q(x) • n(x) dT — J /(x) dfl 

f2e fie 


Cgsb — 


fie 
Nele 

E 

e=l 


de) 

-LSB 


(4.15a) 


(4.15b) 


where c(x) and q(x) are the solutions obtained by solving a given unconstrained LSFEM. In numer¬ 
ical /i-convergence study, we are interested in the following quantities with respect to /i-refinement: 


CMaxAbsLSB — m^ 

U u 0 U u ^ 

EAbsGSB = I^GSbI 


de) I 
^lsbI 


(4.16) 

(4.17) 


Few remarks about the species balance are in order. In writing equation (|4.17p . we have assumed 
that the mesh is conforming, and the test and trial functions belong to C®(n) (i.e., there is inter¬ 
element continuity of the functions). Under the proposed computational framework, we place ex¬ 
plicit (equality) constraints to meet e^gB =0 Ve = 1, • • • ,Nele. By meeting the local species 
balance, the global species balance is trivially met. 


4.4. Discrete equations. Let Kcc denote the stiffness matrix obtained by lower-order finite 
element discretization of the LSFEM terms involving c(x) and w{x). Similarly, we can define the 
stiffness matrices K^q, Kqc, and iUqq- The load vectors are denoted by Vc and Cq, respectively. 
These vectors are obtained from the finite element discretization of the LSFEM terms involving 
u;(x) and p(x). It should be noted that the stiffness matrices Kcc and F^qq are symmetric and 
positive definite. Furthermore, iUqc = K"^^. For more details, see Appendix C. 

The corresponding constrained optimization problem in the discrete setting for the proposed 
locally conservative DMP-preserving LSFEMs can be written as follows: 

1 1 

minimize -(c; KccC) + (c; Kc^q) + -(q; Kq^q) - (c; Vc) - {q; Vq) (4.18a) 

c g Jjrtcdo/S 2 Z 

q g ^nqdofs 

{ AcC -|- AnO = bf , , 

(4.18b) 

Cminl E C E Cmaxl 

where ^‘nqdofs^' denotes the number of degrees-of-freedom for the flux vector, and ^^ncdofs” denotes 
the number of degrees-of-freedom for the concentration. The vector of size ncdofs x 1 with all 
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entries to be unity is denoted as 1. Recall that (•;•) denotes the standard inner-product on the 
Euclidean spaces. The finite element discretization of the local species balance equation gives rise 
to the global LSB matrices and Aq, and the global LSB vector bf. The matrices Ac and Aq 
are of sizes Nele x ncdofs and Nele x nqdofs, respectively. Similar inference can be drawn on 
the sizes of bf, Vc, rq, ITcq, and ITqq. Since iTqc = and the matrices Kcc and iTqq are 
symmetric and positive definite, the constrained optimization problem (|4.18al) - (l4.18bp belongs to 
convex quadratic programming and has a unique global minimizer m The corresponding first- 
order optimality conditions - popularly known as the Karush-Kuhn-Tucker (KKT) conditions - for 
this discrete optimization problem take the following form: 


KccC -\- Kcafl — fc A^ Ac -|- /Jjcin Rmax 

(4.19a) 

■^cq^ T -^qq9 = '*’q “ ^q-^q 

(4.19b) 

AcC -p Aq^ bf 

(4.19c) 

Rmin ^ 0 

(4.19d) 

Rmax — ^ 

(4.19e) 

(c Cjdijil) • /^jnin ^ 

(4.19f) 

(Cmaxl ~ c) • Aljciax ~ 0 

(4.19g) 


where Ac and Aq are the Lagrange multipliers enforcing the LSB equality constraints, which stem 
from equation (I4.19cl) . Hcnin Rmax KKT multipliers enforcing the DMP inequality 

constraints given by Cminl ^ c and c ^ Cmaxl- Note that the non-negative constraint is a subset 
of the DMP inequality constraints. To wit, setting Cmin = 0 and Cmax = +oo will result in explicit 
non-negative constraints on the nodal concentrations. 

Remark 4.1. Note that the continuous problem, in general, cannot always be written as an 
optimization problem. This is eertainly the case with respect to ADR equation m- Moreover, in 
the continuous setting, the non-negative and local species balance constraints are satisfied trivially. 
Therefore, the Lagrange multipliers are zero in the continuous setting (if one can write the problem 
as an optimization problem). The violations of the non-negative and local species balanee constraints 
are only in the discrete setting. This is the reason why one needs to obtain the discrete form before 
one can fix the deficiencies in solving the discrete equations. 

In the next two sections, we illustrate the performance of the proposed computational frame¬ 
work for advection-dominated ADR problems and transport-controlled irreversible fast bimolecular 
reactions. In all the numerical simulations reported in this paper, the constrained optimization 
problem is solved using the MATLAB’s |33| built-in function handler quadprog, which has a robust 
solver based on an interior-point numerical algorithm presented in References [5Tti53|. One can 
alternatively employ the open-source optimization solvers such as COBYLA, SLSQP, L-BFGS-B, or TNG 
from SciPy |54) . The tolerance in the stopping criterion for solving convex quadratic programming 
problems is taken as lOOemach; where Cmach ~ 2.22 x 10“^® is the machine precision for a 64-bit 
machine. 

There are various approaches to numerically solve transient diffusion-type systems. It is desirable 
to have a numerical strategy that converts and utilizes the solvers for steady-state diffusion-type 
equations to solve transient systems. It has been recently shown that the method of horizontal lines 
using the backward Euler time-stepping scheme is one of the viable approaches to respect maximum 
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principles and the non-negative constraint in the discrete setting |55) . The method of horizontal 
lines discretizes the time domain first, and thereby converts the transient ADR equations at each 
time-level into a system of governing equations similar to (|2.1ap - (|2.1cl) . This methodology, thus, 
helps us to use the computational framework provided in Section 01 One can employ a numerical 
procedure similar to Algorithm 1 provided in Reference |55| to advance the numerical solution over 
the time. Numerical results for transient systems are presented in Section [6l 


5. NUMERICAL /i-CONVERGENCE AND BENCHMARK PROBLEMS 


We shall employ a popular problem from the literature, which is commonly used to assess the 
accuracy of numerical formulations for advective-diffusive systems (e.g., see [T4l[56l). The test 
problem is constructed using the method of manufactured solutions. The computational domain 
is a bi-unit square: D = (0,1) x (0,1). The advection velocity vector field is taken as v(x) = e^, 
where is the unit vector along the y-direction. The scalar diffusivity is denoted by D(x). The 
concentration field is taken as follows: 


c{x,y) 


sin(7rx) 

gm2—mi _ 2 




(5.1) 


where the constants mi and m 2 are given in terms of the scalar diffusivity: 


mi 


m2 


1 - Vl -F 47r2D2 
TD 

1 Vl + 47r2D2 
TD 


(5.2a) 

(5.2b) 


We have taken D(x) = 10“^ in our numerical simulations. This choice is arbitrary, and is primarily 
motivated to check whether the proposed framework gives stable, reliable, and accurate numerical 
results for advection-dominated problems. For the chosen value of the diffusivity, the solution (15.ip 
exhibits steep gradients near the boundary of the domain. A pictorial description of the boundary 
value problem is provided by Figure [3j 

Numerical solutions for the concentration and the flux vector are obtained by prescribing Dirich- 
let boundary conditions on all the four sides of the computational domain. These conditions are 
enforced strongly and are given as follows: 


c(x) 


sin(7rx) for y = 0 
0 for x = 0orx = lory=l 


(5.3) 


Using equation (15.ip . one can calculate the corresponding flux vector and volumetric source needed 
for the convergence analysis. 


5.1. Convergence analysis for D{x,y) = 10“^. Herein, we will discuss the performance of 
negatively stabilized streamline diffusion LSFEM with and without LSB constraints. In case of 
unconstrained setting, we also quantify the errors incurred in satisfying LSB and GSB. Numerical 
simulations are performed using a series of hierarchical structured meshes based on three-node 
triangular (T3) and four-node quadrilateral (Q4) elements with XSeed and YSeed ranging from 11 
to 81. Figure [4] provides the typical computational meshes used in the numerical /i-convergence 
analysis. 

The weights for the primitive and negatively stabilized streamline diffusion LSFEMs are taken to 
be of LS Type-1 (i.e., A(x) = I and /3(x) = 1). The element stabilization parameters for negatively 
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stabilized streamline diffusion LSFEM are taken as 6 o = 0.01 and Tq = 0.01. The convergence of 
the proposed computational framework with respect to L 2 -norm and ff^-semi-norm is illustrated 
in Figure O From this figure, one can notice that near optimal convergence rates are achieved for 
the concentration field in both L 2 -norm and Ff^-semi-norm for unconstrained negatively stabilized 
streamline diffusion formulation. For the flux vector, near optimal convergence rate is obtained in 
L 2 -norm but not in iF^-semi-norm. This is because of the steep gradients in the concentration field 
at the boundary y = 1, which is due to the small value for the diffusivity. Enforcing LSB constraints 
considerably improves the iF^-semi-norm convergence rate for the flux vector. However, for the flux 
variables, there is a slight decrease in L 2 -norm convergence rate as compared to the unconstrained 
negatively stabilized streamline diffusion LSFEM. Similar decrease in convergence rates of L 2 -norm 
and FF^-semi-norm for the concentration has been observed. This can be attributed to the fact that 
LSB constraints improve the accuracy of the flux vector inside the boundary layers but has little 
effect away from it. 

Remark 5.1. It should he noted that the convergence rates reported in Figure for the un¬ 
constrained negatively stabilized streamline diffusion LSFEM are in accordance with the mathemat¬ 
ical analysis provided by Kopteva m and Stynes !j26^ . These results are obtained for singularly 
perturbed advection-diffusion equation based on a class of unconstrained streamline diffusion finite 
element formulations. Kopteva m shows that one can get at best first-order convergence inside 
boundary and characteristic layers even on special meshes. 

From Figure [5] one can also conclude that the Q4 element performs better than the T3 element. 
These trends in the convergence rates for different meshes are due to the fact that higher-order 
derivatives (e.g., div[grad[c(x)]) in the stabilization terms for negatively stabilized streamline dif¬ 
fusion LSFEM vanish for T3 element. But these stabilization terms are non-zero for a Q4 element. 
The reason is that the shape functions for a T3 element are affine while that of a Q4 element are 
bilinear. 

Another important aspect of this numerical /i-convergence study is to quantify the errors in¬ 
curred in satisfying LSB and GSB for unconstrained LSFEMs. The contours of the error distribution 
in LSB and the Lagrange multipliers enforcing the LSB constraints are shown in Figure [H It is 
apparent that errors incurred in satisfying LSB are smaller under Q4 meshes than under T3 meshes. 
The decrease in CMaxAbsCSB and CAbsCSB on /i-refinement is shown in Figure [71 From this figure, one 
can notice that the errors in LSB and GSB for a Q4 mesh are lesser than that of a T3 mesh. On 
fi-refinement, the decrease in CMaxAbsLSB and CAbsGSB is slow and not close to machine precision. 

Finally, the computational cost of the unconstrained and constrained LSFEMs for both T3 and 
Q4 meshes are shown in Figures [8] and [9l It is clear that the computational cost associated with a 
Q4 mesh is higher than that of a T3 mesh. This can be again be attributed to the non-vanishing 
stabilization terms (e.g., div[grad[c(x)]) in the negatively stabilized streamline diffusion LSFEM 
for Q4 meshes. For constrained LSFEMs, the maximum additional computational cost (for both 
LSFEMs) did not exceed 15%, which has been tested on a hierarchy of meshes. 

5.2. Thermal boundary layer problem. This benchmark problem has wide practical ap¬ 
plications in the areas of heat and mass transfer. Herein, we shall use this benchmark problem to 
study the performance of unconstrained and constrained LSFEM formulations in capturing steep 
gradients near the boundary for advection-dominated scenarios. Gonsider a rectangular domain 
H = {{x,y) G [0,1] X [0,0.5]} with velocity field v(x,y) = 2yex, where e^, is the unit vector along 
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the x-direction. The volumetric source is assumed to be homogeneous (i.e., f{x,y) = 0), and the 
scalar diffusivity is taken to be D{x,y) = 10“^. The boundary conditions are: 


c{x,y) 


0 for 0 < X < 1 and y = 0 
2y for X = 1 and 0 < y < 0.5 

< 

1 for 0 < X < 1 and y = 1 
1 for X = 0 and 0 < y < 0.5 


(5.4) 


A pictorial description of the boundary value problem is provided in Figure [TOl The weights are 
taken to be that of LS Type-1 (see equations (|4.6ap and ()4.6bl) h The element-level stabilization 
parameters for negatively stabilized streamline diffusion LSFEM are taken to be 5o = 0.01 and 
To = 0.001. Numerical simulations are performed using four-node quadrilateral mesh with XSeed = 
41 and YSeed = 21. The element Peclet number will then be Pe/i = 125. The obtained concentration 
contours are shown in Figure [TTJ It is evident from these figures that numerical solution obtained 
from the primitive LSFEM contains node-to-node spurious oscillations. These oscillations did not 
reduce even after enforcing the LSB and NN constraints. But the negatively stabilized streamline 
diffusion LSEEM is able to capture the steep gradients near the boundary without producing spu¬ 
rious oscillations. The errors incurred in satisfying LSB for unconstrained LSFEM formulations are 
shown in Figure fT^ 


6. TRANSPORT-CONTROLLED BIMOLECULAR CHEMICAL REACTIONS 

In this section, we shall apply the proposed mixed LSFEM-based computational framework 
to study transport-controlled bimolecular chemical reactions. Specifically, we are interested in the 
spatial distribution, plume formation, and chaotic mixing of chemical species at high Peclet numbers. 
To this end, consider the following irreversible bimolecular chemical reaction: 

UA^ + ub B — > nc C (6.1) 

where A, B, and C are the species involved in the chemical reaction; n^, ub, and nc are their 
respective (positive) stoichiometric coefficients. The fate of these chemical species are governed by 
the following coupled advective-diffusive-reactive system: 

-hdiv[vcA-D(x,t)grad[cA]] = fA{x,t) - nArix,t,CA,CB,cc) innx]0,X[ (6.2a) 

dc 

-h div[vcB - D(x,t)grad[cB]] = - nBr{x,t,CA,CB,cc) in Ox]0,X[ (6.2b) 

dc(j 

+ div[vcc - D(x,t)grad[cc]] = /c(x,t) + ncr{x,t,CA,CB,cc) in flx]0,X[ (6.2c) 

Cj(x,t) = cf(x,t) onrfx]0,X[ (6.2d) 

- —v(x,t)ci(x,t) - D(x,t)grad[ci(x,t)]^ • n(x) = hf(x, t) onr^x]0,X[ (6.2e) 

Ci(x, t = 0) = c°(x) in n (6.2f) 

where i = A, B, and C. v(x, t) is the advection velocity vector field, fi{x,t) constitutes the 
non-reactive volumetric source, cf (x, t) is the Dirichlet boundary condition, and h?(x, t) is the 
Neumann boundary condition of the i-th chemical species. r{x,t,CA,CB,cc) is the bimolecular 
chemical reaction rate, which is a non-linear function of the concentrations of the chemical species 
involved in the reaction. c?(x) is the initial condition of i-th chemical species, t S [0,X] denote the 
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time, where X is the total time of interest. The coupled governing equations (|6.2ap - (|6.2el) can be 
converted to a set of uncoupled advection-diffusion equations using the following linear algebraic 
transformation: 


cp := CA + 



(6.3a) 


CG := CB + f^ cc (6.3b) 

\ncj 

As we are interested in fast bimolecular chemical reactions, it is acceptable to assume that the 
chemical species A and B cannot co-exist at any given location x and time t. Hence, c^, cb, and 
cc can be evaluated as follows: 



(6.4a) 


(6.4b) 


(6.4c) 


In Reference |35| . a similar mathematical model has been studied in the context of maximum 
principles and the non-negative constraint. However, the study has neglected the advection, and 
did not address local and global species balance. These aspects are very important and cannot be 
neglected in the numerical simulations of chemically reacting systems. In particular, advection can 
play a predominant role in the study of bioremediation |58| , transverse mixing-controlled chemical 
reactions in hydro-geological media |59| . and contaminant degradation problems |60| . This paper 
precisely addresses such problems in which advection is dominant, and satisfying species balance at 
both local and global levels is extremely important. 

Remark 6.1. Non-linear chemical dynamics is a huge field with various interesting artifacts, 
which include chaos and limit cycles Non-linear reactions will bring many additional 

complications, which need to be addressed systematically. Our approach can handle zeroth-order and 
first-order kinetics, as these two cases do not bring additional challenges. Other reaction kinetic 
models need to be addressed case-by-case. A general treatment of non-linear chemical dynamics is 
not trivial, and is beyond the scope of this paper. 

Herein, we perform numerical simulations for highly spatially varying advection velocity fields 
and time-periodic flows. See Reference |61| for a discussion on time-periodic flows. For such 
problems in 2D, the following quantity is of considerable importance, which is referred as the 
position weighted second moment of the product C concentration: 



/ 


cg(x, t) dD 


where yo is the location of a convenient reference horizontal line. In our numerical simulations, we 
have taken yo to be the y-coordinate of the start of the formation of product C. Since cc{x,t) > 0, 
0^(t) is a non-negative quantity. In subsequent sections, we study the utility of this quantity as a 
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posteriori criterion to assess numerical accuracy. We also analyze the variation of 0^ with respect 
to Pe/i- We also present the numerical results that shed light on the impact of advection on the 
formation of the product C. In all our numerical simulations, we have taken the weights in primitive 
and negatively stabilized streamline diffusion LSFEMs to be that of LS Type-1 (i.e., A(x) = I and 
/3(x) = 1). 

Remark 6.2. In the literature, to study mixing processes due to advection, spectral methods 
pseudo spectral methods m, and model reduction methods m are commonly employed. However, 
such methods are limited to time-periodic flows, periodic initial and boundary conditions, simple ge¬ 
ometries, and homogeneous isotropic diffusivity. Extending these methods to complicated geometries, 
general initial and boundary conditions, complicated advection velocity fields, and heterogeneous 
isotropic and anisotropic difl^usivity is not trivial and may not even be possible. Moreover, these 
methods do not guarantee the satisfaction of non-negativity, DMPs, LSB, and GSB. The proposed 
computational framework is aimed at filling this lacuna. 

6.1. One-dimensional steady-state analysis of product formation in fast reactions. 

Analysis is performed for two different advection velocities: v = 0.25 and v = 1.0. Diffusivity is 
assumed to be 2.5 x 10“^. The stoichiometric coefficients are assumed to be: ua = 2, ub = 1, and 
nc = 1. Numerical simulations are performed for two different cases as described below. 

6.1.1. Case ffl. A pictorial description of the boundary value problem is shown in Figure [TSl 
The objective of this case study is to analyze whether the proposed negatively stabilized streamline 
diffusion LSFEM can produce physically meaningful values for cfix) on coarse meshes. Based on 
the linear algebraic transformation given by equations (I6.3aj) - (l6.3bh . the analytical solution for 
invariants F and G can be written as follows: 


cf{x) 

cg{x) 


/ l — exp(vx/D)\ 

V l-exp(v/D)J 

h f_ 1 - exp(ux/T>) \ 

V k l — exp{v/D) J 


(6.6a) 

(6.6b) 


Using equations (I6.4a|) - (l6.4cjl . one can obtain the analytical solution for product C. 

For the numerical solution, we have taken XSeed = 11. The element stabilization parameters 
for negatively stabilized streamline diffusion LSFEM are taken as 6o = 0.08 and Tq = 0.04 when 
V = 0.25. For v = 1.0, 6o and Tq, are assumed to equal to 0.083 and 0.0121, respectively. The 
analytical and numerical solutions are compared in Figure [TTl As per this figure, the primitive 
LSFEM produces node-to-node oscillations near the boundaries of the domain. Furthermore, its 
numerical solution considerably deviates from the analytical solution in the entire domain. For 
Pe/i = 5 and Pe/j = 20, the negative value for the concentration is as low as —1.25 and —0.47. On 
the other hand, the negatively stabilized streamline diffusion LSFEM is able capture the analytical 
solution profile in the entire domain without producing negative values in the concentration field. 

6.1.2. Case A pictorial description of the boundary value problem is provided in Figure 
M The objective of this case study is to examine whether the proposed LSFEM can capture steep 
gradients in the solution near the boundary. The analytical solution for the invariants F and G can 
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be written as follows: 


cf{x) = 1 - 


1 — ex.p{vx / D) 


cg(x)= 


1 — exp{v/D) 
1 — exp{vx / D) 


(6.7a) 


(6.7b) 


1 — exp{v/D) 

Figure [15] compares the obtained the numerical solution with the analytical solution. The negatively 
stabilized streamline diffusion LSFEM is able to accurately capture the steep gradients near the 
boundary. 


6.2. Steady-state plume formation from boundary in a reaction tank. A pictorial 
description of the boundary value problem is provided in Figure [T6| The computational domain is 
a rectangle with = 2 and Ly = 1. Dirichlet boundary conditions with = Cg = 1 are specified 
on the left side of the domain. Elsewhere, c?(x) is taken to be zero for all the chemical species 
involved in the bimolecular reaction. The non-reactive volumetric source is assumed to be zero in 
the entire domain for all the chemical species. The stoichiometric coefficients are taken as = 1, 
ns = 1 and nc = 1- The advection velocity field is defined through the following multi-mode 
stream function |35| : 


where x= {x,y), {pi,P 2 ,P 3 ) = (4,5,10), (gi, 92 , 93 ) = (1,5,10), and (^ 1 ,^ 2 , A 3 ) = (0.08,0.02,0.01). 
The corresponding components of the advection velocity can be written as follows: 

3 


r ^ 1 , fPkTTx 7r\ /9fc7ry\ 

' + - 2 ) [-Li) 

''.w = +757 = (— - ^™ (—) 


k=l 


\ J 


(6.9a) 

(6.9b) 


It is easy to check that div[v(x)] = 0. The contours of the stream function and the corresponding 
advection velocity vector field are shown in Figure [16] Numerical simulations are performed using 
the following two different types of diffusivities: 

• Type #1: D{x) = 10“^ 

• Type ^2'. D(x) = RDqR^, where R and Dq are given as follows: 


R = 


Do(x) = ujo 


cos(0) — sin(0) \ 
sin(0) cos(0) ) 

f yl + (^ 2 x 1 

\ -(1 - 0J2)x^y^ 


-(1 - uj2)x*y* \ 
^22/* +xl ) 


(6.10a) 

(6.10b) 


where = x + oji and 9 * = y + oji. The parameters 0, cjq, Wi, and U 2 are equal to 
1.0, 10“^, and 10“^. Correspondingly, the eigenvalues of D(x) are ujq{x‘ 1 + yl) and 
0 JQUJ 2 {xt + 9 *). The contrast/anisotropic ratio of the media (which is the ratio of maximum 
to minimum eigenvalue) is as high as 10^. 

Herein, we employed a structured mesh based on Q4 elements. Numerical simulations are performed 
with varying mesh sizes and polynomial orders (p = 1,2,3) to demonstrate the pros and cons 
of various unconstrained and constrained LSFEMs. The stabilization parameters are taken as 
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So = To = 10“^ and S 2 = T 2 = 10“^. The contours of the concentration of the product C are shown 
in Figures [T7H201 for both the primitive and negatively stabilized streamline diffusion LSFEMs. 
The white patches in the figures denote the regions in which the non-negative constraint has been 
violated. The variation of 0^ with respect to XSeed and Pbl are shown in Figures [2Tfl22l From 
these figures, the following inferences can be drawn: 

(i) ft is clear that both low-order and higher-order polynomials violate the non-negative constraint 
and DMPs under unconstrained formulations. Moreover, mesh refinement and polynomial 
refinement do not seem to reduce the amount of violated region for DMP constraints. 

(ii) The proposed framework based on p = 1 is able to satisfy all the desired properties, and is 
able to predict physically meaningful values for the concentration and the flux. 

(iii) The primitive LSFEM and the unconstrained negatively stabilized streamline diffusion LSFEM 
give unphysical values for the position weighted second moment of the product C (i.e., 0^). 
On the other hand, the proposed computational framework is able to accurately describe the 
variation of 0^ with respect to mesh refinement. In addition, the numerical values for 0^ 
reaches a plateau on /i-refinement, which indicates convergence. However, this is not observed 
with the unconstrained primitive and negatively stabilized streamline diffusion LSFEMs. 

Finally, it should be emphasized that placing explicit non-negative constraints on the nodal concen¬ 
trations does not ensure non-negativity of the concentration in the entire computational domain. 
This is due to the fact that higher-order shape functions change their sign within an element |66) . 


6.3. Transient analysis of non-chaotic and chaotic vortex stirred mixing in a reac¬ 
tion tank. Figure [231 provides a pictorial description of the problem with appropriate initial and 
boundary conditions. The computational domain is a square with = Ly = 1. For all chemical 
species, zero flux boundary condition is prescribed on the entire boundary. The non-reactive vol¬ 
umetric source is zero in the entire domain for all the chemical species A, B, and C. Reactant A 
is placed at the center of vortices, which are positioned at (0.25,0.75) and (0.75,0.25). The width 
of the square slug A is equal to 0.25. Within this slug, CA(x,t = 0) = 8. Elsewhere, the initial 
condition for A is equal to zero. Correspondingly, the initial condition for reactant B is zero in 
these two square regions centered at (0.25,0.75) and (0.75,0.25). Elsewhere, CB(x,t = 0) = 1.5. 
The stoichiometric coefficients are taken as ua = I, = I, and nc = 1- The total time of interest 
is taken as X = 5. We assume scalar diffusivity to be D = 10~^. The stabilization parameters 
are taken as So = To = 10“^ and Si = Ti = 10““^. For advection velocity, we employ the following 
vortex-based flow field: 

• Type : Non-chaotic vortex-based advection velocity field 


v(x) = cos(27ry)ea; -t- cos(27rx)ey 
Type ^2\ Chaotic vortex-based advection velocity field 

cos(27ry) -|- Vo sin(27ry) if vT <t< [u + ^)T 


V3;(x,t) = 


cos(27ry) 


if (i/+i)r<t < (i/+l)r 


( 6 . 11 ) 


( 6 . 12 ) 


Vy(x,t) = 


cos(27ra;) 

cos(27rx) -|- Vo sin(27rx) 
27 


if i/T < t < [ly + T 
if {iy + l)T <t < {u + l)T 


(6.13) 


where u = 0,1,2, ■■ ■ mm- T denotes the period of the motion of the flow field. Vq is an a-priorly 
chosen chaotic flow perturbation parameter. Herein, we choose T = 0.8 and Vq = 1-0 m- 

Figures [21H23 provide the concentration profiles of unconstrained and constrained negatively 
stabilized streamline diffusion LSFEM with NN constraints. Herein, XSeed = YSeed = 121. Nu¬ 
merical simulations are performed for various different time steps. These are equal to 0.0001, 0.001, 
0.01, and 0.1. Figure [2^ shows the concentration profile of the product C and Figure [25] shows 
the values of Cc at y = 0.5 at the first time-step. Analysis is performed using the unconstrained 
weighted negatively stabilized streamline diffusion LSFEM. From these figures, it is clear that un¬ 
physical negative values for cc are obtained even for small time steps. Furthermore, these violations 
are significant and not close to machine precision for both small and large time-steps. Non-negative 
constraints have to be enforced in order to get meaningful values for cc- 

Figures [26] and [27| show the concentration profiles of cc for both non-chaotic and chaotic vortex 
flow fields at various time levels. For non-chaotic advection, product C is initially formed away from 
the vortex field. As time progresses, it slowly gets accumulated in the closed streamlines of the two 
vortices. Regions of higher concentration are located at the center of vortices. From Figure [26l it is 
evident that cc contour is symmetric along the line y = x. This is because the non-chaotic vortex- 
based advection velocity vector field is symmetric along this line. However, this is not the case with 
chaotic vortex flow field. Qualitatively, there is no symmetry associated with the concentration 
field. This is because of the time-periodic sinusoidal terms given by equations (I6.12h and (I6.13h . 
They provide chaotic features for the chosen value of period of motion T. An interesting feature 
observed is that mixing of chemical species is enhanced in chaotic flow as compared to non-chaotic 
flow. This is because cc is not equal to zero in the non-chaotic flow for late times. 

Finally, from these figures it is evident that existing numerical formulations do not provide 
accurate information on the fate of reactants and products for all times. On the other hand, the 
proposed methodology predicts results accurately for both early and late times. 


6.4. Transient analysis of species mixing in cellular flows. A pictorial description of the 
initial boundary value problem is provided in Figure[28] The prescribed diffusive/total flux for each 
chemical species is taken to be equal to zero. The initial condition is such that c^(x) = 1 in the 
bottom half of the domain and vanishes elsewhere while Cg(x) = 1 only in the upper half. The 
non-reactive part of the volumetric source is equal to zero for all the species. It should be noted 
that the concentration of the product C should be between 0 and 1 because = 0. 

For numerical simulations, we have taken Lx = 1, Ly = 0.5, XSeed = 61 and YSeed = 241. 
The stoichiometric coefficients are taken as ha = 1, ns = 1, and nc = 1- The time-step is taken 
as At = 0.1. The total time of interest is taken as X = 5. The scalar diffusivity is taken as 
D = 5 X 10“^. The stabilization parameters are taken as do = Tq = 10“^ and 6i = Ti = 10“^. The 
advection velocity vector field for the cellular flow is given by |68| : 


V X = — sm 


27rx 


L 


Cell 


COS 


27ry 

f^Cell 


ex + cos 


27rx 


L 


Cell 


sm 


/ 2TTy 
V -^Cell 


(6.14) 


where Lceii is the cell length of a pair of vortices. The velocity field given by equation (I6.14h has 
a set of symmetrical vortices, and the neighboring vortices rotate in opposite directions. It is well- 
known that the advection velocity field given by equation (I6.14|) causes numerical difficulties (if the 
underlying numerical scheme is not properly designed). This is because the advection velocity is 
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strongly non-uniform [681 Section 8], (which happens to be in our case). That is, 

X At > 1 (6.15) 

The main objective of this test problem is show that the proposed formulation is robust and can ana¬ 
lyze velocity fields that are strongly non-uniform without causing numerical instabilities/oscillations. 

Analysis is performed for a series of hierarchical cell lengths. That is Tceii equal to 0.5, 0.25, 
0.125, 0.0625, and 0.03125. In all the cases, as the input data and the position of the cellular 
vortices is symmetric about the line y = 0.25, it may he expected that formation of product C will 
be symmetric along this line. Additional information on the symmetry of the formation of product 
C can be inferred based on the cell length of adjacent pair of vortices, which rotate in opposite 
directions. This is apparent from the numerical results presented for product C provided in Figures 

MM 

Figure [29] shows the concentration profiles of the product C under the unconstrained and con¬ 
strained negatively stabilized streamline diffusion for Lceii = 0.5. From this figure, it is evident that 
the unconstrained LSFEM violates both the non-negative and maximum constraints. In addition, 
both undershoots and overshoots are observed. On the other hand, the proposed computational 
framework with LSB and DMP constraints provides physically meaningful profiles for the concen¬ 
tration of the product C. 

For Tceii = 0.5, immediately after time t = 0, we observe wing-like concentration profiles. This 
is because diffusion controls the species mixing rather than advection across the adjacent cells along 
the line y = 0.25 (note that Vy{x,y = 0.25) = 0). However, once the species A and B enter the 
closed streamlines where advection dominates, mixing happens at much faster rate. Furthermore, 
product C spreads in time along the array of counter-rotating vortices. After a considerable time 
(t ~ 5), we observe that formation of product C is symmetric along the lines x = 0.5 and y = 0.25. 
In addition, product C accumulates near the region where ||v(x)|| is close to zero (which happens 
to be at the center of vortices and hyperbolic points). This happens to be in a good agreement with 
the inferences drawn from the numerical simulations performed on cellular flows |61[l68j . 

Figure [29] and [3T] show that the separatrices connecting the hyperbolic points inhibit long range 
transport of chemical species from one cell to another. By decreasing Lceii, species mixing can be 
enhanced. This is evident from Figure [3TJ Qualitatively, the numerical results presented here agree 
with the analysis presented by Neufeld and Garcia m, which tells us that mixing of chemical 
species is fast within a cell but the transport of reactants/products between the cells is controlled 
by diffusion only. In order to enhance the efficient species mixing in different regions in these type 
of flows, Lceii has to be as smaller. To conclude, we would like to emphasize that the numerical 
solution based on the proposed methodology does not exhibit numerical instabilities and is able to 
capture the essential features even when the advection velocity is strongly non-uniform. 

7. SUMMARY AND CONCLUDING REMARKS 

We presented a robust computational framework for (steady-state and transient) advection- 
diffusion-reaction equations that satisfies the non-negative constraint, maximum principles, local 
species balance, and global species balance. The framework can handle general computational 
grids, anisotropic diffusivity, highly heterogeneous velocity fields, and provides physically meaningful 
numerical solutions without node-to-node spurious oscillations even on coarse computational meshes. 
The main contributions of the paper can be summarized as follows: 

29 


max 


dvx 

dx 


dy 




(Cl) We constructed and proved a continuous maximum principle that includes both Dirichlet and 
Neumann boundary conditions. It also takes into account the inflow and outflow Neumann 
boundary conditions in establishing the maximum principle. 

(C2) We described in detail the shortcomings of several plausible numerical approaches to satisfy 
the maximum principle, the non-negative constraint, and species balance. 

(C3) We proposed a locally conservative DMP-preserving computational framework and constructed 
element stabilization parameters that are valid for a general reaction coefficient, advection 
velocity, and diffusivity. The framework has been carefully constructed using the least-squares 
finite element method (LSFEM). It is also shown that a naive implementation of LSFEM will 
not meet the desired properties. 

(C4) The discrete problem under the proposed framework is well-posed, and it can be shown that 
a unique solution exists. 

(C5) We performed numerical convergence studies on the computational framework. We also sys¬ 
tematically analyzed and documented the performance of the proposed framework with various 
benchmark problems and realistic examples. 

(C6) We obtained numerically a scaling law for a transport-controlled bimolecular reaction. 

(C7) In chemically reactive systems, it is important to predict the fate of reactants and products 
during the early times. We have shown that the existing formulations may not provide accurate 
information for such scenarios. On the other hand, using numerical experiments, we have 
shown that the proposed framework predicts accurate results for both early and late times. 


The salient features and performance of the proposed computational framework can be summa¬ 
rized as follows: 

(51) The rate of decrease of errors in LSB and GSB for unconstrained negatively stabilized stream¬ 
line diffusion LSFEM with /i-refinement is slow and is about 0{h). Furthermore, this numer¬ 
ical formulation violates various discrete principles and the non-negative constraint for both 
isotropic and anisotropic diffusivities. On the other hand, the proposed non-negative compu¬ 
tational framework is able to satisfy LSB and GSB up to machine precision on an arbitrary 
computational mesh. 

(52) The proposed computational framework with NN and LSB constraints eliminates the spurious 
node-to-node oscillations and provides physically meaningful values for concentration. Fur¬ 
thermore, it is able to furnish reasonable answers with various time-steps and at various time 
levels even on coarse computational grids. 

(53) It has been shown that existing formulations may fail to give acceptable results for non¬ 
negative statistical quantities such as 0^, which is defined in Section 16.21 However, the 
proposed methodology always provides non-negative values for 0^. The quantity 0^ can be 
used as a posteriori criteria to assess accuracy of numerical solutions for complex initial and 
boundary value problems for which non-negativity and species conservation are important. 

(54) Due to the aforementioned desired properties, our proposed computational framework can be 
an ideal candidate to numerically obtain scaling laws for complicated problems with non-trivial 
initial and boundary conditions. Therefore, the proposed framework will be vital for predictive 
simulations in groundwater modeling, reactive transport, environmental fluid mechanics, and 
modeling of degradation of materials. 
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A possible future research work is to implement and analyze the performance of the proposed 
numerical methodology in a parallel environment. A related research is to design tailored iterative 
solvers and associated pre-conditioners for our proposed numerical methodology. 


8. APPENDIX A: Element-level discretization of stabilization terms 

The weighted negatively stabilized streamline diffusion least-squares formulation requires the 
evaluation of div[grad[c(x)]] and grad[grad[c(x)]] terms at the element-level. Since these terms are 
not typical, we present a compact way of discretizing these terms under the finite element method. 
It should be noted that one need not evaluate these terms for lowest-order simplicial elements (i.e., 
three-node triangular element and four-node tetrahedron element), as these terms will be identically 
zero. However, div[grad[c(x)]] and grad[grad[c(x)]] will be non-zero for high-order simplicial and 
non-simplicial elements (i.e., four-node quadrilateral element and eight-node brick element). 

8.1. Fourth-order tensors. Let R and S be two second-order tensors. A fourth-order tensor 
product R Kl S is defined as follows: 

(RKS)T = RTS'^ (8.1) 

for any second-order tensor T. The fourth-order identity tensor can then be written as: 

I = IKI (8.2) 

where I is the second-order identity tensor. The fourth-order transposer T and symmetrizer S 
tensors are defined as follows: 

TA = A'^ (8.3a) 

SA = ^ (A + A^) (8.3b) 

where A is a general second-order tensor. 


8.2. Properties of Kronecker products relevant to finite element discretization. Let 

A and B be matrices of sizes n x m and p x q, respectively. Let Uij and bij be, respectively, the 
ijth components of A and B. The Kronecker product of these two matrices is an np x mq matrix 
that is defined as follows: 


AQB : = 


auB . 

• • ^ImB 

^nlB . ■ 



Another operator that will be useful is the vec[«] operator, which is defined as follows: 


(8.4) 


On 


vec[A] := 




O'nl 


(^nm 


(8.5) 
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The following standard properties can be established for mat[«] and vec[«] operators: 


vec [A + B] = vec [A] -|- vec[P] 

(8.6) 

(A QB){CQD) = (AC 0 BD) 

(8.7) 

vec[ACP] = (P^ 0 A) vec[C] 

(8.8) 


For general properties of Kronecker products, see the book by Graham |69) . However, this book 
does not contain many of the results presented below, which are useful for the current study. 

For an effective computer implementation of LSFEM-based formulations, we shall represent 
four-dimensional and three-dimensional arrays as two-dimensional matrices. To this end, consider 
a four-dimensional array P of size m x n x p x q. The matrix representation of P is denoted by 
mat[«], and is defined as follows: 



Piiii • 

Pllpl 

Plll2 

■ Plllq 

mat[P] := 

Pmnll 

Pmnpl 

Pmnl2 

P 

• mnpq 


The following properties can be established for the matrix representation of fourth-order tensors: 


vec [PX] = mat [P] vec [X] 
mat[A M B] = B Q A 

mat[§] = ^ (7 0 / -|- mat[T]) 

where the matrix representation of T takes the following form: 


mat[T] : = 


J0[l,0,0,...,0]„ 

J0[0,l,0,...,0]„ 


J0[O,O,O,...,1]„ 


mnxmn 


( 8 . 10 ) 

( 8 . 11 ) 

( 8 . 12 ) 


(8.13) 


In the above equation, I is an identity matrix of size m x m. It can be shown that for a matrix Z 
of size m X n, vec[Z"'"] = mat[T]vec[Z]. 

For a three-dimensional arrays, there are two useful matrix representations, which will be de¬ 
noted as mati[«] and mat 2 [«]. Consider a three-dimensional array Q of size m x n x p. The 
corresponding matrix representations of Q are defined as follows: 


mati[Q] := 


Qlll 

Qlnl ■ ■ ■ 

• • • Qllp 

Qlnp 






(8.14) 

Qvrill 

Qmnl ■ ■ ■ 

• • • Qmlp 

Qmnp 
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mat2[Q] := 


<5iii 

• • • Qmll 

Qiip 

• • • Qmlp 

Qi21 

• • • Qm21 

Ql2p 

• • • Qm2p 

_ Qlnp 

• • • Qranp _ 


(8.15) 


The following properties can be established for the matrix representations of three-dimensional 
arrays: 


vec[( 51 ^] = mati[Q]vec['F] (8.16) 

Yec[Qz\ = mat 2 [Q]vec[ 2 :] (8-17) 


8.3. Finite element discretization of div[grad[c(x)]] and grad[grad[c(x)]] terms. Let ^ 
denote the position vector in the reference finite element. The row vector containing the shape 
functions is denoted by N , which is a function of The derivatives and Hessian of N with respect 
to ^ are, respectively, denoted as DN and DDN . That is, in indicial notation we have: 

{DN\^ = ^ (8.18a) 

where Ni and are, respectively, the T’^-component of the vectors N and The concentration 
field c(x) and total flux vector q(x) are interpolated as follows: 

c(x) = (^(®)) (8.19a) 

q(x) = (^(a:)) (8.19b) 

_^rp _^rp 

where c and q denote the matrix containing nodal concentration and total flux vector. In 
negatively stabilized streamline diffusion LSFEM, div[grad[c(x)]] and grad[grad[c(x)]] terms arise 
from the following expansions: 

div[L)(x)grad[c(x)]] = grad[T)(x)] • grad[c(x)] -|- Il(x) div[grad[c(x)]] (8.20a) 

div[D(x)grad[c(x)]] = div[D(x)] • grad[c(x)] -|- D(x) • grad[grad[c(x)]] (8.20b) 

Based on the regularity of diffusivity, grad[Z)(x)] and div[D(x)] can be calculated in multiple ways. If 
the diffusivity is continuously differentiable, then grad[Il(x)] and div[D(x)] can be directly evaluated 
analytically. This will considerably reduce the computational cost in the evaluation of local stiffness 
matrices. If D{x) is not differentiable (but is square integrable), then grad[Z)(x)] and div[D(x)] 
can be evaluated as follows: 

grad[T»(x)] = {DN)J-^ (8.21a) 

div[D(x)] = mati[T)]vec[mat[T](T)A7') J“^] (8.21b) 
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where J is the Jacobian matrix, and D is the nodal values for the diffusivity, whose size can be 
inferred based on the context (whether the diffusivity is either -D(x) or D(x)). Using equations 
(18.181) and (|8.19p . the Laplacian and Hessian of c(x) can be calculated as follows: 

div[grad[c(x)]] = “ {DN)J-^x^^ mati[i:>i:>Af]vec[J-^ vec[c^] (8.22) 

vec[grad[grad[c(x)]]] = ^ mat 2 AT] [l Q I — xj'^{DN)'^) ^ vecfc"*"] (8.23) 

where x is the matrix containing nodal coordinates. 


9. APPENDIX B: Coercivity, error estimates, and stabilization parameters 


Herein, we shall establish coercivity and error estimates for the homogeneous weighted LSFEMs 
(i.e., cP(x) = 0 on did). Based on this mathematical analysis, we obtain the stabilization parameters 
that are used in this paper for the weighted negatively stabilized streamline diffusion LSFEM. To 
this end, let denote the standard Sobolev space for a given non-negative integer m |28| . 

The associated standard inner product and norm are denoted by (•;«),^ and || • \\m, respectively. 
On the function space Q, the inner product (•;«)(jiv and norm || • ||div are defined as follows: 

(p; q)div := (p; q)o + (div[p]; div[q])o Vp, q £ Q (9.1a) 

IIpIIL := IIpIIo + l|div[p]||o Vp G Q (9.1b) 

The Poincare-Friedrichs inequality takes the following form |45| : there exists a constant Cpf > 0 
such that we have 


■alio < C'p/||grad['u]||o Vm G 


(9.2) 


Consider the classical weighted primitive least-squares functional S^prim ((c, q) , /) given by equa¬ 
tion (|4.5p with c(x) = 0 on 912. If (c, q) G C X Q is an exact solution of the equations (|2.1ap - (|2.1cp . 
then (c, q) must be a unique zero minimizer of Sprim ((c, q) , /) on C x Q. Hence, for any e G M, we 
have 


A. 

de ' 


5'Prim((c, q)+ e(^n,p) ,/) =0 V (w, p) G W X Q 


e=0 


which is identical to the following: 


(9.3) 


'Bprim ((C, q) ; {W, p)) = Tprim {{w, p)) 


V (tc, p) G >V X Q 


(9.4) 


where iBprim(*;*) and Tprim(*) are the corresponding bilinear and linear forms for the weighted 
primitive least-squares functional Jprim- It should be noted that 


)Bprim ((■;«, p) ; (tf^,p)) = SPrim ((t(^,p) ,/ = 0) V (w, p) G W X Q (9.5) 

Equation (19.5p is used to prove coercivity and boundedness estimates for the bilinear form fBprim- 
Now consider the finite element discretization of the equation (j9.4p . Let Ch © C, JV/i C W, and 
Qh © Q be the hnite element function spaces spanned by piecewise polynomials of degree less than 
or equal to r over 12/^. It should be noted that r is an integer and r > 1. Then, the discrete weighted 
primitive LSFEM can be written as follows: Find (c/i,q/i) G Ch x Qh such that 

53prim ((ch, q/i) ; {Wh, Ph)) = -Cprim {{wh, Ph)) V {wh, Ph) S Wh X Qh 
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(9.6) 











where (c/i,q/i) is the finite element solution with respect to the chosen basis functions spanning the 
finite element space Ch x Q^. Similar inference holds for 5NgStb((c, q),/), ^NgStb; and -CNgStb- 
We assume that is quasi-uniform mM- That is, there exists a constant C > 0 independent 
of h such that h < for all Qg G Additionally, we assume that the following inverse 

inequality holds on these quasi-uniform meshes. There exists a constant C > 0 independent of h 
such that 


^ div[grad[c,j]] < ||grad[c/,]||g ^Ch^Ch 


(9.7a) 


D • grad [grad [cft]] 


< 


tr[D]tr[grad[grad[c/j]]] ^ = tr[D] div[grad[cfe]] ^ Vc/i G Ch (9.7b) 


where tr[«] is the trace of a matrix. In proposing equation ()9.7bl) . we assumed that the Hessian of 
Ch, grad [grad [cft]], is positive semi-definite. 

All the results presented here are applicable for a general anisotropic diffusivity tensor, advection 
velocity vector field, and linear reaction coefficient. One can obtain simplified results for isotropy 
by taking D(x) = L)(x)I, where I is an identity tensor. 


Theorem 9.1 (Coercivity for weighted primitive LSFEM). There exist constants Cprimi > 
0 and Cpr\ra 2 > 0 independent o/D and h such that for all {wh,Ph) G x Qh-' 


T?Prim (('R’/i! P/i) )/ — 0) > (7pi.inil7niinAjjj;^||grad[u;/i] IIq 


T?Prim ((R’/i) Ph) , f — 0) ^ C*Prim27min'^ 
where the positive constant 7min is: 

7min := min 


2 

min 




It-. I|2 

I11 div 


1 + 


(9.8a) 

(9.8b) 

(9.9) 


1, min [/3(x)] , min [Amin,A(x)] 
x&n xefi 

where Aniin,A(x) is the minimum eigenvalue o/A(x) at a given point x G id. 

Proof. Consider the weighted primitive least-squares functional (14.Sp with / = 0. Equation 
(j9.9p implies: 


25pr 


„2 > 


Ph - Wh'v + Dgrad[u;,i] - g.gTa.<i[wh] 


o,n 


-k 


awh + div[p/j] - fiWh 


o,Q 


+ 2/r (ph - WhV + Dgrad[u;/,]; grad[u;/i])o + div[ph]; u;/i)o 

- T^\\wh\\l,n - fi^\\gra.d[wh]\\ln 


(9.10) 


where /r is a positive constant, which will be determined later. Using PoincarAFriedrichs inequality 
and Green’s formulae, equation (|9.10p can be written as: 


25pr 


X 


> fi (2Amin - (l + Cp/)) IlgradKlIlg^tj 


We obtain equation (I9.8ap by choosing 


h = 


1 + ^P/ 


(9.11) 


(9.12) 
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There exist two non-negative constants and Ca (for instance, Cy, = [||v||^] and Ca = 

xGfr 

Ta^l) such that 

x&n 


l|w^/^lll = llw^ftllo + IlgradKlIlo < 




7^ ■ 

I mm mm 


-5^Pr 


IPhllo < 2||p/j - rnftv Dgradfin^] ^ + 2|| - Wh'v + Dgrad[u;fe] 


0,C2 


< 1 + 


2C^C^f(l + C^f) x2 >4^^ 

1 +- + + 


, M „2 „ „2 / + 

||div[p/,]||o < 2||au;h-gdiv[p/,]||g^^-g2||au;fe||Q^^ < IH-- 

\ ^min 


45^Pr 


(9.13a) 


7„ 


(9.13b) 

(9.13c) 


It is easy to check that inequalities (|9.13al) - (l9.13cl) imply inequality (I9.8bp . □ 

Theorem 9.2 (Coercivity and boundedness estimate for NSSD LSFEM). Given that equa¬ 
tions (|9.7ap - (I9.7bp hold. If for each Hg S we take 


4 {ndfXl^^ + CCl^5a,y,h? + 

__ 

32 (l + Clf) + CClfd^^h? + C5^h^) 


(9.14a) 

(9.14b) 


then for all (wh, Ph) £ VV/i x Qh there exist two constants C^NgStbo > 0 and CNgStb 4 > 0 independent 
o/D and h such that we have: 

Coercivity 


((»*. P.) . / = 0) > 


+ E 


c7L.tii„'‘L 


I grad[u;/; 


llo,ne 


32(1 + Cl^) (nd2AL. + CC^,fSa.h^ + C6j,hA 


(9.15) 


Boundedness estimate 

CNgStbi||R’h||i + C'NgStb2||p/i||div + C^NgStbsllv • grad[u;h]||o ^ S^NgStb {{wh, Ph ), / = 0) 

< C'Ngstb4 (||r^/i||i + llp/illdiv + l|v« grad[u;,j]||o) 
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(9.16) 


























where the constant 7min is given by equation M . The eonstants ^NgStbi; ^NgStb 2 ? and 

C*NgStb 2 given as follows: 


Sav 

<5d 


CNgStbl 

C'NgStb2 


C'NgStbS 


(a + div[v])^ 


xGH •- 

[||div[D]|p] 

xGH 

„2 \2 

'min 


C*NgStb07min-^r 


_^NgStb07min^min*^avD_ 

(l + + A^ax) '^ovD + ^«vDA^jjj/i^ + (5 q,viA^;^/i^ 

C*NgStb07min'^min^ 

^avD 


(9.17a) 

(9.17b) 

(9.17c) 

(9.17d) 

(9.17e) 


The constants 5avi and SavT> in the above equations are defined as follows: 


(5„vi = max 


xGf2 
'max 


(grad [a] • v + adiv[v])^ 


<^q:vD — "^max d^vh -|- djjh 


(9.18a) 

(9.18b) 


Proof. The boundedness estimate is a direct consequence of the triangle inequality. Herein, 
we shall proceed to show the validity of coercivity estimates, specifically, equation (|9.15p and the 
left hand side of (j9.16p . Let /r > 0 be a constant, which will be determined later. Using equation 
(19.9p and (14.131) with / = 0, we have 


2 

0,Oe 


2gNgStb ^ I p/j - rc/jv + Dgrad[r(;ft,] - (5neV(div[rc/iV - Dgrad[rc/i]]) -//gradfrch] 

I 2 

+ ^ \ awh + dw[ph\ + 6n^dW[awh^] - fiwh - f?\\wh\\l^^ - ir^\\grad[wh\\\l^Q 




+ ^ awh + div[wh'v-DgTcad[wh]] ^ ^ 2fi {awh + div[ph + dn^awh^r]-Wh)Q^Q 

+ '^ (.Ph - Wh^ + Dgrad[wfe] - (J^eV (div[u;/jV - Dgradfu;/*]]); grad[wfe])Q (9.19) 


37 

























Using Theorem 19.11 equation (I9.14ap - ()9.14bl) . Cauchy-Schwartz inequality, Poincare-Friedrichs in¬ 
equality, Green’s formulae, and following inequalities 


2Tn^ {(a + div[v]) Wh]^ • grad[u;fe])o 
-2rne ((a div[v]) Wh; D • grad[grad[u;h]])o 
-2rne (v • grad[wfe]; D • grad[grad[u;h]])o 
2rn, ((a div[v]) Wh', div[D] • grad[u;/,])o 
-2 to^ (div[D] • grad[u;fe]; D • grad[grad[u;h]])o 
2rne (v • grad[wfe]; div[D] • grad[wfe])o 


> TnJ\ (a-hdiv[v])u;,,||g_f2^ 


+ rnJ|v*gradK]||^f^^ 

(9.20a) 

> rnell (a + div[v])u;/,||g_f2^ 


+ T-HellD • grad[grad[u;/j]]||g f^^ 

(9.20b) 

> rnJ|v*gradK]||^ t^^ 


+ TneljD • grad[grad[u;h]]||o,ne 

(9.20c) 

> TnJ\ (a + div[v])u;/,||g_f2^ 


+ rnJ|div[D] • grad[u;fe]||§^t 2 ^ 

(9.20d) 

> rnJ|div[D] .gradKill 


+ TneljD • grad[grad[u;/j]]||o,n. 

(9.20e) 

> rnJ|v*gradK]||^f^^ 


+ rnJ|div[D] • grad[u;fe]||g_f 2 ^ 

(9.20f) 


we have the following inequality: 


m, ate/idiv[u;,jv - Dgrad[u;ft,]] 


E 


> — 


0,Ce 


16 I 1 -|- Cpj 


• grad[wh]\\l^^ 


16(1 + C^f) (nd^Xl^ + + Cd-oh^) 

Similarly, using the following equality: 


(9.21) 


2/rfe, {diY[awhv];wh)o^Q^ = -2/ifee {awh] v • grad[u;,j])o (div[av]u;/,; (9.22) 


in combination with the following inequalities 


-2//fe^ ((a -h div[v]) Wh] v • grad[u;/,])o 
2/ife^ (div[D] • gradPu;/,]; V • grad[u;/i])o 
2/ife, (D • grad[grad[u;fe]]; V • grad[u;/i])o^Q^ 


> 2/i(5f2ell (« + div[v])u;/,||g o 


_l_ 11 ^ ^ gradfrc- 


0,f1e 


> 2/i(5f2el|div[D] •grad[u;ft]||§^f2^ 


+ • gradfrc 


0,Qe 


> 2/i(5f2j|D • grad [grad [re, 
+ ^^l|v*gradK]||go 


0,Oe 


(9.23a) 


(9.23b) 
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(9.23c) 















and choosing fj, = , equation (|9.19p reduces to the following: 


-'pf 


25^NgStb ^ 3A 


2 

min 


+ E 




4 (l + 8(1 + (nd^A^^ + CCl^S^^K^ + C6^h?) 


2 

0,Qe 


(9.24) 


+ ^ arcft + div[t(;/iv - Dgradftc/i]] 

From equations (19.211) and (|9.25ah . we get the desired result given by equation (19.151) . The second 
part of the proof is similar to Theorem 19.11 These exist a constant Cq-vD > 0 (for instance, 


CavD = max 


nd^,C,CC^f 


such that 


Ajnax + ~\~ CSjyh ^ C^avD^QvD 


||gradK]||g < 


|v* gradf-WhJIlo < 


32(1 + C^f)dNgSth 
llyT A2. 

/mm mm 

32C'ovDdavD(l + C'p/)C'^T?NgStb 


Using Cauchy-Schwartz inequality on ||v • grad[u;ft,]||o and ()9.25bp gives 


2 ||2|| ir i||2 32 C'v( 1 + C'pj)5^NgStb 


|v*grad[u;fe]||^ < ||v||g||grad[u;fe]||^ < 
Now, consider the terms ||tCh||i and ||ph|ldiv 

lo + llgrad[u;/,]||g < 


lly^ • A^ 

I min min 


2 II or ll|2 ^ 32(1 + C'pj)^5NgStb 


llo^ ■ A^ ■ 

I min min 

Ph\\l<2 ^ Ph - WhV + Dgradlwh] - (div[r(;/iV - Dgrad[r(;ft,]]) 


(9.25a) 

(9.25b) 

(9.25c) 

(9.26) 

(9.27a) 




2 

O,0e 


+ 2 ^ - Wh^r Dgrad[r(;fe] - (div[rc/iV - Dgrad[rc/i]]) 


n^efih 


2 

o,ne 


(9.27b) 


||div[p,,]||2 < 2 ^ awh + div[pfe] + (5oediv[atc/iv] 




2 

O,0e 


2 

0,Oe 


(9.27c) 


+ 2^1 awh + (5oediv[ar(;/iv] 

OeSOh 

Using (I9.25ap - (j9.26p and repeated use of triangle inequality on ()9.27bp and (j9.27cp gives the bound¬ 
edness estimate (I9.16F □ 

Theorem 9.3 (Error estimate for proposed LSFEM). Given that equations (I2.1al) - (|2.1cl) have 
a sufficiently smooth solution (c, q) G (C x Q) n (i7'’'*'^(Al)) . Then the finite element solution 
(c/uQ/i) of the unconstrained weighted negatively stabilized streamline diffusion LSFEM satisfies the 
following error estimate: 


V C'NgStblllc — Ch\\l + V C'NgStb2||q — qh,||div + V ^NgStbsllv • grad[c — C/j]||o 

E C^gStb/i (ll^llr-l-l T l|q||r-|-l) 

where CNgStb > 0 is a constant independent o/D and h. 
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(9.28) 























































Proof. Let c/ G Ch and q/ G Qh be the standard finite element interpolants of c and q, 
respectively. From the interpolation theory |45| . we have 

||c — c/||i < C'/i'’||c||r+i (9.29a) 

||q - q/lldiv < C’/i’’||q||r-+i (9.29b) 

for some positive constant C independent of D and h. The error (c — c/i, q — q/i) satisfies the 
following orthogonality property: 

®NgStb ((c/i - c,q/i - q); (rt;/i,ph)) = 0 V (wfe, p/*) G W/i x (9.30) 

Cauchy-Schwartz inequality implies: 

®NgStb - Ci,(lh - q/) ; {ch - C/, q/j - q/)) < ^B^gstb ((c - c/, q - q/) ; (c - c/, q - q/)) (9.31) 

From Theorem 19.21 and interpolation estimates (I9.29ap - (j9.29bh . one can obtain the desired error 

estimate (19.281) . □ 


From the above mathematical analysis, it is evident that the element-dependent stabilization 
parameters rog < 0 and < 0 can be taken as: 


Sfie = - 


= 


^o^minhQ^ 


+ <^1 niax (a -I- div[v])^ /i ^-|-^2 [||div[D] |p] 

xgO 1 f xGO 

2 u2 


T 


^max + n max 

xGO 


(a-|-div[v])^ /i^-|-r 2 m^ [||div[D]||^] h 


(9.32a) 

(9.32b) 


where 5o, hi, 62 , Tq, ti, and T 2 are non-negative constants. 


Remark 9.4. The mathematical analysis provided by Hsieh and Yang can be obtained as 
a special case of the mathematical analysis presented above. Specifically, take a = 0, D(x) to be 
homogeneous and isotropic, and v(x) to be solenoidal and constant. 


10. APPENDIX C: Finite element stiffness matrices and load vectors 


For weighted primitive LSFEM the terms and Cq are constructed from the 

local stiffness matrices and load vectors IT^q, -f^qqi 'I'c: and r® through the standard hnite 

element assembly process. The expressions for these element stiffness matrices and element load 
vectors in terms of shape functions and their derivatives are explicitly defined as follows: 




-1 

fie 


N^Ndne + J dQe - 

fie f1, 

-IfT 


J (DA1J~^) DA^vAl dfle 


- J AT^v^A^D (DNJ-i)" dfle + J {DNJ-^)'DA^'D{DNJ-^)^ dflg 

fie fie 

2 TvrT ; 


+ 


1 -|- Sign[v • n] 


{vfiyN'^NdTl 


( 10 . 1 ) 


fieer? 
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= 

cq 


/ {/?"<») 


AT'"'" (vec 


£>Arj-^)^j)^ dQe - I AT^v^A^ (AT 0 J) d!^e 


+ / (DAT J-i) DA2 (AT 0 J) 

Oe 

1 + Sign[v • n] 


Cesr'? 


(v • n) (AT 0 I) dF^ 


= J ^2 (^vec [(£)Arj"i)^j j (^vec [(£>Arj"i)'^'J dOe + J {N^ Q l) {N Q I) dF!e 
6 0 
+ j {N^ Q I) {N Q I) dTl (10.3) 

Oeer'j 

Correspondingly, the expressions for the element load vectors in terms of shape functions and their 
derivatives are explicitly defined as follows: 

1 + Sign[v • n] 


( 10 . 2 ) 


= 


He 

< = / {/3V) 


dOe - 


vec 


Oesr? 

-ilT 


(v • n) AT^gP dr^ 


(10.4) 

(10.5) 


(Z>AT'J-^)"J dfle + J {N^Ql)nqPdri 

Sle SleSr'? 

It should be noted that these terms are obtained from the bilinear and linear forms of the weighted 
primitive least-squares functional ^primj which are given as follows: 

53prim ((c, q) ; {w, p)) = (w; ^‘^a^c) + (w; (v • A^v) c) - (gradfrc]; (DA^v) c) 

— (w;v • A^Dgradfc]) -|- (grad[r(;]; (DA^D) grad[c]) 

+ (u;;/32adiv[q]) -|- (div[p];/il^ac) — (u;;v« A^q) — (p;A2vc) 

+ (grad[r(;]; DA2 q) + (p; A^D grad[c]) -|- (div[p];/32div[q]) 


+ (p; A^q) + 


w: 


1 + Sign[v • fi] 


(v • n)^ c 


r? 


- w; 


1 + Sign[v • n] 


(v • fi) q • n I + (p • n; q • fi) 


r<j 


1 -|- Sign[v • fi] 


p*n; 

£pnm((w,p)) = (w;/3^af) + (div[p];/32/) - ( w; 

+ (p«n;gP; 


r? 

T'? 

1 -|- Sign[v • fi] 


(v • fi) c 


( 10 . 6 ) 


(v • fi) gf 


r<2 


)r, (10.7) 

Similarly, one can derive the stiffness matrices and load vectors for weighted negatively stabilized 
streamline diffusion LSFEM S^NgStb- Foi’ sake of saving space, herein we shall not explicitly dehne 
them as the bilinear and linear forms of S^NgStb have more than fifty terms (from which Kcq, 

Vc, and Vq are derived). 
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0.015 



Figure 2. Academic problem: This figure compares the numerical solutions under the stan¬ 
dard single-field Galerkin formulation and the normal equations approach with the exact 
solution. The normal equations approach does not eliminate node-to-node spurious oscilla¬ 
tions. 


y 



Figure 3. Numerical h-convergence study: A pictorial description of the two-dimensional 
boundary value problem used in the numerical convergence analysis. Dirichlet boundary 
conditions are prescribed on the entire boundary. 
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log(error) , log(error) 



(a) Mesh using T3 elements (b) Mesh using Q4 elements 


Figure 4. Numerical /i-convergence study: This figure shows the typical computational 
meshes used in the numerical convergence analysis. The meshes shown in this figure have 21 
nodes along each side of the computational domain (i.e., XSeed = YSeed = 21). A series of 
hierarchical computational meshes are employed in the study with 11 x 11, 21 x 21, 41 x 41 
and 81 x 81 nodes. 
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(a) Concentration: No constraints 



(b) Concentration: LSB constraints 




Figure 5. Numerical h-convergence study: This figure shows the convergence rates for the 
concentration and flux vector in L 2 -norm and iJ^-semi-norm with and without LSB con¬ 
straints. Convergence studies are performed using T3- and Q4-based meshes under the 
negatively stabilized streamline diffusion LSFEM. It is evident that the Q4 element slightly 
outperforms the T3 element in terms of rates of convergence. 
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(a) T3 mesh: Error in LSB 


(b) T3 mesh: Lagrange multiplier en¬ 
forcing LSB 
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(c) Q4 mesh: Error in LSB 


(d) Q4 mesh: 
forcing LSB 


Lagrange multiplier en- 


Figure 6 . Numerical h-convergence study: The top and bottom left figures show the con¬ 
tours of error incurred in satisfying LSB for unconstrained negatively stabilized streamline 
diffusion LSFEM. The right set of figures show the contours of Lagrange multiplier enforcing 
LSB constraint using the proposed LSFEM. Note that the Lagrange multipliers enforcing 
the LSB constraint can have negative value as opposed to KKT multipliers. Numerical 
simulations are performed based on three-node triangular mesh and four-node quadrilateral 
mesh with 81 nodes on each side of the domain. In essence, the LSB errors and Lagrange 
multipliers enforcing LSB based on a Q4 mesh is lesser than a T3 mesh. 
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Figure 7. Numerical /i-convergence study: These figures show the decrease of CMaxAbsLSB 
and CAbsGSB with respect to XSeed for a series of hierarchical three-node triangular and four- 
node quadrilateral meshes. (See equations (I4.15I) - (I4.17I) for the definitions of CMaxAbsLSB and 
CAbsGSB-) Numerical simulations are performed using the unconstrained primitive and nega¬ 
tively stabilized streamline diffusion LSFEMs. For XSeed = 81, CMaxAbsLSB and CAbsGSB are 
in 0(10“®). In addition, the decrease in LSB and GSB errors with respect to /i-refinement 
is slow, and the values are not close to the machine precision. 




Figure 8 . Numerical h-convergence study: This figure shows the CPU time (in seconds) 
of the proposed computational framework for unconstrained primitive and unconstrained 
negatively stabilized streamline diffusion LSFEMs. For Q4 mesh, as div[grad[c]] 0, the 
computational cost is higher than that of the T3 mesh. 
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Computational overhead (%) 



(a) T3 mesh 



Figure 9. Numerical /i-convergence study: This figure shows the computational overhead 
incurred in satisfying LSB as compared to that of the corresponding unconstrained formu¬ 
lations. Analysis is performed for primitive and negatively stabilized streamline diffusion 
LSFEMs. For XSeed = 11, we obtained negative value for the computational overhead. This 
is because the interior point convex algorithm used in MATLAB optimization solver |33| 
pre-processes the constrained convex quadratic programming problem simplifies the given 
LSB constraints by removing redundancies. Hence, for very low number of unknowns, 
the computational cost associated with interior point convex algorithm is much faster 
than the LU solver for the unconstrained optimization problem. 

(0,0.5) c=l (1,0.5) 



Figure 10. Thermal boundary layer problem: This figure shows a pictorial description of the 
boundary value problem. Dirichlet boundary conditions are prescribed on all four sides of 
the computational domain. We have taken c(x) = 1 at x = (0,0). 
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(a) Primitive (No constraints) 



0.800 

0.700 

0.600 

0.500 

0.400 

0.300 

0.200 

0.100 


0.00 


0.00 


(b) Negatively stabilized streamline diffusion (LSB 
and NN constraints) 


Figure 11. Thermal boundary layer problem: This figure shows the contours of concentration 
obtained for both unconstrained and constrained LSFEMs based on Q4 finite element mesh. 
The proposed LSFEM-based framework with NN and LSB constraints is able to eliminate 
spurious oscillations near the boundaries y = 0 and x = 1. This is not the case with the 
primitive LSFEM. 
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(b) Negatively stabilized streamline diffusion 


Figure 12. Thermal boundary problem: This figure shows the contours of the error incurred 
in satisfying LSB for various unconstrained LSFEM formulations using Q4 meshes. One 
can notice that the error is more dominant in the interior of the domain under the primitive 
LSFEM, whereas the error is dominant at the boundary x = 1 under the negatively stabilized 
streamline diffusion formulation. 


cP (x = 0) = 1 
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Ca(x = 1) = 0 

/c = 0 
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cP (x = 1) = 0 


Figure 13. ID irreversible bimolecular fast reaction problem: A pictorial description of the 
boundary value problem. For Case ^1: /b(x) = 1 and Cg(x = 1) = 0, and for Case ^2: 
fB{x.) = 0 and c^(x = 1) = 1. 


50 












cc(x) ca(x) 



Figure 14. ID irreversible bimolecular fast reaction problem (Case #1): This figure com¬ 
pares the concentration profile of the reactants and the product for various element Peclet 
numbers using unconstrained primitive and unconstrained negatively stabilized streamline 
diffusion LSFEMs to that of the analytical solution. The primitive LSFEM considerably 
deviated from the analytical solution. Moreover, it violated the non-negative and maximum 
constraints. On the other hand, the negatively stabilized streamline diffusion LSFEM is 
able to capture the analytical solution exactly in the entire domain even at high element 
Peclet numbers. 
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Figure 15. ID irreversible bimolecular fast reaction problem (Case #2): This figure compares 
the concentration profile of the chemical species A, B, and C for various element Peclet 
numbers using unconstrained primitive and unconstrained negatively stabilized streamline 
diffusion LSFEMs to that of the analytical solution. The negatively stabilized streamline 
diffusion LSFEM is able to capture the features near the boundary layer with considerable 
accuracy even on coarse meshes. 
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cP(x) = 0 



(a) Problem description 



(b) Stream function and advection velocity vector field 

Figure 16. Plume development from boundary in a reaction tank: The top figure provides a 
pictorial description of the boundary value problem. The bottom hgure shows the contours 
of the stream function corresponding to the advection velocity vector held. 
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Figure 17. Plume development from boundary in a reaction tank (Type #1): This figure 
shows the concentration profiles of the product C based on unconstrained primitive LSFEM. 
The white region shows the area in which concentration is negative. Both the lower-order 
and higher-order finite elements violate the non-negative constraint. The negative values 
are in the range 0(10“^) to 0(10“'^), which are not close to the machine precision Cmach = 
0 ( 10 - 1 ®). 



0 0 0.5 1.0 1.6 2 0 


(c) XSeed = YSeed = 251 (LSB and NN constraints) 

Figure 18. Plume development from boundary in a reaction tank (Type #1): This figure 
shows the concentration profiles of the product C based on unconstrained and constrained 
negatively stabilized streamline diffusion LSFEM. The white region shows the area in which 
concentration is negative. Considerable part of the domain violated the non-negative con¬ 
straint. The proposed framework with NN and LSB constraints is able to capture the plume 
formation on a coarse mesh for a highly heterogeneous advection velocity vector field. 
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(a) p = 1, XSeed = YSeed = 101 


(b) p = 1, XSeed = YSeed = 501 



Figure 19. Plume development from boundary in a reaction tank (Type #2): This figure 
shows the concentration profiles of the product C based on unconstrained primitive LSFEM. 
The white region indicates the area in which the obtained concentration is negative. The 
negative values are in the range 0(10“^) to 0(10“®). Both h-refinement and p-refinement 
could not eliminate the violation in the non-negative constraint for this problem, which has 
highly heterogeneous anisotropic diffusivity. 


55 











(c) XSeed = YSeed = 251 (LSB and NN constraints) 


Figure 20. Plume development from boundary in a reaction tank (Type #2): This figure 
shows the concentration profiles of the product C based on unconstrained and constrained 
negatively stabilized streamline diffusion LSFEM. Compared to Figure [1^1 the proposed 
methodology with NN and LSB constraints is able to accurately describe the plume for¬ 
mation of the product C even for highly heterogeneous anisotropic diffusivities and highly 
spatially varying velocity fields. 



(a) No constraints 



(b) With LSB and DMP constraints 


Figure 21. Plume development from boundary in a reaction tank (Type #1): This figure shows 
the variation 0^ with mesh refinement under the weighted negatively stabilized streamline 
diffusion LSFEM. It should be noted that 0^ is a non-negative quantity. However, the 
unconstrained negatively stabilized streamline diffusion LSFEM gives negative values for 
0p. Mesh refinement did not alleviate this problem. On the other hand, the proposed 
framework not only gives non-negative values for 0^ but flattens upon mesh refinement, 
which indicates convergence. 
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Figure 22. Plume development from boundary in a reaction tank (Type #1): This figure 
shows the variation log(0^) with respect to v/Pcl for isotropic diffusivity under the weighted 
negatively stabilized streamline diffusion LSFEM with LSB and DMP constraints. Herein, 
analysis is performed using XSeed = YSeed = 201. Through numerical simulations, we 
observed that log(0^) oc VFei. 
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/lP(x,t) 


hi (x, i) = 0 




(a) Problem description 



(b) Stream function and associated advection velocity field 



(c) Reactant A: Initial condition 


(d) Reactant B-. Initial condition 


Figure 23. Vortex-stirred mixing in a reaction tank: The top-left figure provides a pictorial 
description of the initial boundary value problem. The top-right figure shows the contours 
of the stream function corresponding to the advection velocity field given by (16.111) . The 
bottom figures show the initial conditions for the reactants A and B such that = 

0)) = (cs(x,t = 0)) = 1. 
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(c) Product C at t = 0.01 


Figure 24. Non-chaotic vortex-stirred mixing in a reaction tank: This figure shows the con¬ 
centration profiles of the product C after the first time-step using the unconstrained weighted 
negatively stabilized streamline diffusion LSFEM. We have taken XSeed = YSeed = 121. 
If constraints are not enforced, one gets unphysical negative values for the concentration of 
product C. This will be particularly true in the early times of a numerical simulation. 
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(c) cc{y ~ 0.5, t — 0.01) with At = 0.01 




Figure 25. Non-chaotic vortex-stir red mixing in a reaction tank: This figure shows the con¬ 
centration profiles of the product C y = 0.5 after the first time-step using the uncon¬ 
strained weighted negatively stabilized streamline diffusion LSFEM. The violations of the 
non-negative constraint are significant, and are present for various choices of the time-step. 
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(a) Product C at t = 0.1 (No constraints) 


(b) Product C at t = 0.1 (NN constraints) 



(c) Product C &t t = 0.5 (No constraints) (d) Product C at t = 0.5 (NN constraints) 



(e) Product C t = 1.0 (No constraints) (f) Product C at t = 1.0 (NN constraints) 



(g) Product C at t = 5.0 (No constraints) (h) Product Catt = 5.0 (NN constraints) 


Figure 26. Non-chaotic vortex-stirred mixing in a reaction tank: This figure shows the con¬ 
centration profiles of the product C at various time levels using the weighted negatively 
stabilized streamline diffusion LSFEM with and without constraints. The time-step is taken 
as At = 0.1. Herein, XSeed = YSeed = 121. As t increases, the product C should accumu¬ 
late near the center of the two vortices. The proposed computational framework is able to 
accurately capture such features, and the obtj^ined solutions are physical at all times. 











(a) Product C at t = 0.1 



(b) Product C at t = 0.5 





(f) Product (7 at t = 3.0 
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(h) Product C at t = 5.0 


Figure 27. Chaotic vortex-stirred mixing in a reaction tank: This figure shows the concen¬ 
tration profiles of the product C at various time levels using the proposed method with 
non-negative constraints. The time-step is taken as At = 0.1. Herein, XSeed = YSeed 
= 121. An interesting feature observed is that the mixing of cc is enhanced in the entire 
domain when compared to non-chaotic advection (as the minimum value for cc is greater 
than zero). 62 
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(a) Problem description 
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(b) Stream function and advection velocity 

Figure 28. Transport-controlled mixing in cellular flows: A pictorial description of the initial 
boundary value problem and associated advection velocity field for the cellular flow when 
Aceii = 0.5. 
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(a) Product C at t = 0.1 (No constraints) (b) Product C at t = 0.1 (LSB and DMP con¬ 

straints) 



(c) Product (7 at t = 0.5 (No constraints) 
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(d) Product C at t = 0.5 (LSB and DMP con¬ 
straints) 



(e) Product C at t = 1.0 (No constraints) 


(f) Product C at t = 1.0 (LSB and DMP con¬ 
straints) 



(g) Product C &% t = 5.0 (No constraints) (h) Product C at t = 5.0 (LSB and DMP con¬ 

straints) 


Figure 29. Transport-controlled mixing in cellular flows: This figure shows the concentration 
profiles of the product C at various time levels using the unconstrained and constrained 
weighted negatively stabilized streamline diffusion LSFEM when Lceii = 0.5. The proposed 
computational framework is able to produce physically meaningful solution (i.e., satisfies 
the non-negative constraint, maximum principle, and local species balance) for the product 
(7 in a transient cellular flow. The time-step At taken for the numerical simulation is equal 
to 0.1. 
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(a) Product C for Lcell = 0-25 



Figure 30. Transport-controlled mixing in cellular flows (with hierarchical cell lengths): This 
figure shows the concentration prohles of the product C at t = 1 using the proposed for¬ 
mulation with LSB and DMP constraints. Analysis is performed for a series of hierarchical 
Lceii- The time-step At is taken to be equal to 0.1. 
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(b) (cc> vs Lcell 


Figure 31. Transport-controlled mixing in cellular flows (with hierarchical cell lengths): This 
figure shows the average concentration of the product C at various times and different 
cell lengths. Analysis is performed using the proposed formulation with LSB and DMP 
constraints. The time-step At is taken to be equal to 0.1. The main inference from this 
numerical simulation is that species mixing happens faster as Lceii decreases. 
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